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SUMMARY 

The  purpose  of  this  research  is  to  model  a synchronous  generator, 
using  statistical  estimation  theory  to  determine  the  parameters  of  the 
model  from  experimental  data.  An  ideal  generator  model  is  represented 
in  state-space  formulation,  and  the  method  of  quasilinearization  is  used 
to  develop  an  optimal  parameter  estimation  algorithm,  imp  1 emeu tod  on  a 
small  digital  computer. 

The  experiment  to  produce  data  for  the  parameter  estimator  is 
designed  with  the  aid  of  a computer  simulation  of  the  experimental 
setup.  The  simulation  includes  a model  of  a typical  generator  with  a 
switched  load.  The  objectives  of  the  design  of  the  experiment  are  to 
produce  sufficient  data  to  permit  the  estimator  to  work  well,  while 
recognizing  economic  constraints  on  equipment  and  computer  time. 

The  experiment  thus  designed,  is  Implemented  on  a real  synchron- 
ous generator  with  data  recorded  and  digitized  off-line.  The  data  is 
then  stored  on  a magnetic  disk  cartridge  in  a form  suitable  for  use  by 
the  digital  computer.  The  parameter  estimates  enable  the  generator  to 
be  modeled.  The  adequacy  of  the  model  is  validated  by  predicting  res- 
ponse of  the  generator  to  a larger  change  in  load.  The  predicted 
response  matches  the  actual  response  within  the  variance  o t the 
measurement  noise. 

The  standard  methods  of  machine  parameter  estimation  have  been 
formulated  primarily  to  provide  data  for  utility  system  stability 
studies.  The  models  employed  for  this  purpose  are  extracted  from  rather 
simplified  measurements,  requiring  simplifying  assumptions,  and  do  not 
account  for  measurement  noise  and  inaccuracies.  This  report  presents  two 
methods  of  estimating  the  generator  model  parameters  from  sw i t ched- load 
test  data  corrupted  by  considerable  measurement  noise. 


iv 


TABLE  OF  CONTENTS 


SECTION 


INTRODUCTION 


Motivation 

Review  of  Past  Approaches 
Definition  of  the  Problem 


DEVELOPMENT  OF  THE  ESTIMATOR 


2.1  Introduction 

2.2  Cenerator  Model 

2.3  Weighted  Least-Squares  Estimation 

2. A Statistical  Estimation  Theory 

2.4.1  Bayesian  Approach 

2.4.2  Maximum  A Posteriori  Probability  Estimator 

2.4.3  Lower  Bound  on  the  Error  Covariance 

2.5  Summary 


DESIGN  OF  THE  EXPERIMENT 


Introduction 

Equipment  Constraints 

Simulation  Using  Nominal  Parameters 


IMPLEMENTATION  OF  THE  EXPERIMENT 

4.1  Introduction 

4.2  Data  Processing  Methods 

4.3  Instrumentation 

4.4  Result  of  the  Experiment 


CONCLUSIONS. 


5.1  Summary  of  Results 

5.2  Significance  of  Results 

5.3  Recommendations  for  Future  Research 

ESTIMATION  OF  IEEE  STANDARD  PARAMETERS  . . 


6.1  Introduction 

6.2  Parameter  Estimation 

6.2.1  Weighted  Least-Squares  (WLS)  Method 

6.2.2  Maximum  Likelihood  (ML)  Method 

6.3  Application  to  Synchronous  Machines 

6.3.1  Machine  Model  at  Sudden  Short-Circuit 

6.3.2  Machine  Parameter  Estimator 


v 


TABl.E  OF  CONTENTS  (Cont'd) 


SECTION  PAGE 

6.4  Simulation  Results 

6.4.1  WLS  Estimation 

6.4.2  Output  Noise  Effect 

6.4.3  Ml,  Estimation 

6.4.4  Input  Noise  Effect 

6.4.5  Effect  on  Short-Circuit  Resistance 

6.5  Conclusions 

APPENDIX  A LOWER  BOUND  ON  ERROR  COVARIANCE  9 7 

APPENDIX  B NOMINAL  PARAMETERS  FOR  SIMULATION  101 

APPENDIX  C PARAMETER  IDENTIFICATION  PROGRAM  107 

APPENDIX  D GENERATOR  SIMULATION  PROGRAM  123 

APPENDIX  F,  STANDARD  PARAMETER  ESTIMATION  PROGRAM 131 


REFERENCES 


139 


- - 


LIST  OF  ILLUSTRATIONS 

FIGURE  PAGE 

la  Schematic  Cross-Section  of  a Synchronous 

Generator  8 

lb  D-Q  Circuit  Model  of  Synchronous  Generator  Due  to 

Park 10 

2 Schematic  Diagram  of  Load  Bank 29 

3a  Parameter  Error  for  Sudden  Short-Circuit,  Test,  One, 

Parameters  1-5 33 

3b  Parameter  Error  for  Sudden  Short-Circuit,  Test  One, 

Parameters  6-11 34 

3c  Armature  and  Field  Current  for  Sudden  Short-Circuit  . . 35 

4a  Parameter  Error  for  Switched  Resistive  Load,  Test  Two, 

Parameters  1-5  37 

4b  Parameter  Error  for  Switched  Resistive  Load,  Test  Two, 

Parameters  6-11 38 

4c  Armature  and  Field  Current  for  Switched  Resistive  Load, 

Test  Two 39 

5a  Parameter  Error  Versus  Iteration  for  Test  Three, 

Parameters  1-5 40 

5b  Parameter  Error  Versus  Iteration  for  Test  Three, 

Parameters  6-11 41 

5c  Armature  and  Field  Current  for  Test  Three 42 

6a  Parameter  Error  Versus  Iteration  for  Test  Four, 

Parameters  1-5 44 

6b  Parameter  Error  Versus  Iteration  for  Test  Four, 

Parameters  6-11 45 

6c  Armature  and  Field  Current  for  Test  Four 46 

7 Schematic  Diagram  of  Instrumentation  for  Channel  i . . 48 

8 Four  Channels  of  Data  Digitized  and  Written  to  Disk 

Serially 49 


vii 


LIST  OF  ILLUSTRATIONS  (Cont'd) 


FIGURE  PAGE 

9 Reduction  of  Two  Channels  of  Digitized  Data  to  One 

Serial  Record  51 

10a  Parameter  Estimates  Versus  Iteration  Number  from 

Experimental  Data  (Test  One),  Parameters  1,  2,  3,  5, 

6 , and  7 56 

10b  Parameter  Estimates  Versus  Iteration  Number  from 

Experimental  Data  (Test  One),  Parameters  4,  8,  9,  10, 

and  11 57 

11a  Phase  A and  Field  Currents  for  Test  One,  Iteration 

Zero 58 

lib  Phase  A and  Field  Currents  for  Test  One,  Iteration 

13 59 

12  Phase  A Armature  and  Field  Currents  60 

13  Block  Diagram  for  Machine  Parameter  Estimation  ....  79 

14  Oscillogram  of  Short-Circuit  Currents  80 

15  Ratio  of  WLS  Estimate  to  True  Value  for  X’,' 83 

d 

16  Ratio  of  WLS  Estimate  to  True  Value  for  X? 83 

17  Ratios  of  WLS  Estimate  to  True  Value  for  Initial 

Estimates  at  Random 84 

18  Ratio  of  WLS  Estimate  to  True  Voltage  for  Random 

Initial  Estimates:  Output  Noise  Effect  . 87 

19  Ratio  of  ML  Estimate  to  True  Value  for  Random  Initial 

Estimates 8^ 

20  Ratio  of  WLS  Estimate  to  True  Value  for  Random  Initial 

Estimates:  Input  and  Output  Noise  Effect  92 

21  Ratio  of  ML  Estimate  to  True  Value  for  Random  Initial 

Estimates:  Input  Noise  Effect  93 


vlii 


LIST  OF  TABLES 


TABLE  PAGE 

1 Nominal  Parameters  In  Per  Unit 32 

2 List  of  Instruments  Used S3 

3 Parameter  Error  Weights  5S 

4 WLS  Estimation 82 

5 Output  Noise  Effect  on  WLS  Estimate  86 

6 ML  Estimation 88 

7 Input  Noise  Effect  on  WLS  and  ML  Estimates 91 

8 Effect  of  Short-Circuit  Resistance  on  WLS  Estimates  . 94 


ix 


SECTION  I 


INTRODUCTION 


1.1  MOTIVATION 

In  developing  electric  generation  equipment  to  meet  changing  require- 
ments such  as  those  posed  by  military  applications  where  the  machine  is 
used  to  power  a variety  non-standard  loads,  it  is  essential  that  the  machine 
and  load  be  compatible.  It  is  essential,  therefore,  that  techniques  to 
analyze  various  generator  configurations  in  specific  loading  environments 
be  available.  Such  techniques  can  be  implemented  using  digital  computer 
simulation  if  the  parameters  of  mathematical  models  of  the  machines  are 
known.  Therefore,  a problem  of  great  practical  importance  is  the  deter- 
mination of  the  parameters  of  models  of  synchronous  generators. 

The  parameters  of  an  electrical  machine  can  be  calculated,  at  least 
in  principle,  if  the  internal  geometry  and  material  properties  are  known. 
However,  such  detailed  information  is  often  not  available  in  a form  amen- 
able to  system  analysis.  In  this  case,  one  must  subject  the  machine  to 
tests,  recording  data  from  measurements  at  the  terminals.  This  data  is 
invariably  corrupted  with  noise  during  recording  and  processing.  Therefore, 
the  basic  problem  of  determining  machine  parameters  from  tests  is  to  de- 
velop an  estimation  algorithm  to  extract  the  parameters  from  noisy  data 
and  to  design  an  experiment  that  produces  sufficient  data  for  this 
algorithm  to  work  accurately. 

1.2  REVIEW  OF  PAST  APPROACHES 

A coupled  circuit  model  of  a three-phase  synchronous  generator 
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derived  by  Park  [1]  is  in  general  used  by  power  system  analysts.  This 
model  is  obtained  from  a linear  lumped-parameter  model  by  transforming  the 
armature  quantities  onto  a two-axis  coordinate  frame  that  rotates  with  the 
rotor.  Rotor  circuits  consisting  of  a field  winding  and  two  damper  wind- 
ings are  invariant  under  the  transformation.  The  result  of  Park's  trans- 
formation is  a set  of  stationary  differential  equations  in  the  two-axis 
coordinate  frame,  related  to  the  terminal  voltages  and  currents  by  a set 
of  time-varying  measurement  equations.  The  parameters  of  this  model  are 
the  inductances  and  resistances  of  the  two-axis  circuits. 

Since  these  parameters  include  mutual  inductance  between  rotor  and 

i 

stator  circuits  and  inductances  of  rotor  circuits  which  are  not  accessible, 
a simplified  parameter  set  is  often  used  [2,3].  The  simplified  parameters 
can  be  determined  from  results  of  relatively  simple  tests  [4].  An  approxi- 
mate analysis  is  carried  out  by  assuming  that  a transient  in  the  armature 
current  initially  affects  the  damper  windings  and  then  later  affects  the 
f i •=•  1 d winding.  Also,  the  time  constants  of  the  dampers  are  assumed  to  be 
much  shorter  than  the  time  constant  of  the  field.  These  assumptions  are 
generally  accepted  but  tend  to  force  the  characteristics  of  a standard 
machine  on  all  machines  [2]. 

Several  attempts  to  overcome  the  drawbacks  of  this  simplified 
analysis  have  been  published  recently.  Canay  [5]  added  a parameter  to 
account  for  a different  coupling  between  the  various  rotor  and  armature 
circuits.  The  results  predicted  rotor  quantities  with  greater  accuracy. 

Yu  and  Moussa  [b[  then  described  several  approaches  for  determining 
Canay 's  reactance  from  tests.  Transfer  functions  for  synchronous 
generators  have  been  determined  from  sinusoidal  perturbations  about  an 


1 


operating  point  which  are  introduced  by  a fast-response  static  exciter  [7] 
and  from  low-voltage  measurements  with  the  rotor  at  standstill  [8].  The 
first  method  requires  a source  of  field  current  (an  exciter)  which  has  a 
very  fast  response  time  while  the  second  requires  a variable  frequency 
source  of  rather  large  current.  In  these  frequency  response  tests.  Bode 
plot  construction  techniques  were  used.  This  method  requires  judgment  of 
the  analyst  in  locating  the  breakpoint,  since  it  is  essentially  a graphical 
method.  No  treatment  of  errors  was  given.  Stanton  [9]  presented  a sta- 
tistical treatment  of  estimating  a transform  function  that  was  an  empirical 
relation  between  rotor  speed  and  electric  power  output.  The  result  is  not 
directly  applicable  to  determining  the  parameters  of  a physical  model.  Lee 
and  Tan  [10]  recently  implemented  a least-squares  algorithm  to  estimate  the 
parameters  of  the  simplified  analysis  plus  Canay's  reactance.  This  approach 
assumed  that  the  generator  was  subjected  to  a sudden  three-phase  short 
circuit  test.  The  resulting  estimator  was  tested  on  simulated  data  without 
the  effects  of  measurement  noise. 

1.3  DEFINITION  OF  THE  PROBLEM 

The  problem  considered  in  this  research  is  the  estimation  of  the 
inductance  and  resistance  parameters  of  Park's  model  of  a synchronous 
generator.  No  intermediate  parameter  set  requiring  additional  assumptions 
about  the  machine  characteristics  is  imposed.  Data  is  taken  from  the 
terminals  of  a synchronous  generator  under  a resistive  load.  A transient 
is  induced  by  a sudden  switching  of  the  load  resistance  to  a smaller  value. 

A statistical  formulation  of  a parameter  estimation  algorithm  is  studied 
to  enable  an  essessment  of  the  effect  of  measurement  noise  upon  the 
experiment.  This  approach  overcomes  many  of  the  objections  to  previous 


approaches  by  applying  estimation  theory  in  a simple  experiment  to  deter- 
mine the  parameters  of  a synchronous  generator  model  directly. 

A state-space  approach  is  taken  by  casting  the  model  in  the  form 

" f (x.y.u.t)  (1) 

z = h(x,y,t,t)  + w , (2) 

where  x is  the  state  vector,  u is  the  input  vector,  y is  the  parameter 
vector,  z is  the  measurement  vector  and  w is  an  additive  noise  term.  The 
problem  is  to  estimate  y based  on  measurements  z. 

The  first  step  is  to  derive  a parameter  estimator.  This  is  an  al- 
gorithm implemented  on  a digital  computer  for  recursively  computing  esti- 
mates of  the  model  parameters.  A weighted  least-squares  approach  is  taken 
to  minimize  an  error  criterion  that  is  the  weighted  sum  of  the  square  of 
the  error.  The  error  is  defined  as  the  difference  between  the  observed 
output  and  the  output  computed  from  the  model  using  the  current  parameter 
estimates.  An  optimal  estimator  is  derived  from  this  formulation  by 
choosing  the  weighting  matrix  as  the  inverse  of  the  noise  covariance 
matrix. 

The  second  step  is  to  design  an  experiment  which  produces  sufficient 
data  to  permit  the  estimator  to  work  effectively.  The  success  of  the 
parameter  estimator  is  directly  dependent  upon  the  experimental  conditions. 
In  particular,  the  load  and  excitations  ns  well  as  the  data  sampling  rate 
must  be  adequately  chosen.  To  approach  this  problem,  an  experiment  was 
designed  with  the  aid  of  computer  simulations.  That  is,  a simulation  of  a 
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generator  with  nominal  parameters  is  used  to  produce  data  for  the  parameter 


estimator.  This  enabled  the  experimental  conditions  to  be  designed  with 
maximum  flexibility. 

Finally,  an  experiment  was  Implemented  on  an  actual  generator. 

Data  recorded  on  an  FM  Instrumentation  tape  recorder,  was  digitized  and 
put  into  a form  suitable  for  use  by  the  estimation  algorithm.  The  parameter 
estimator  was  then  used  to  estimate  the  generator  parameters.  These  esti- 
mates were  used  to  predict  the  response  of  the  generator  to  a larger  change 
in  load  resistance.  This  final  step  validated  the  results  and  enabled  the 
overall  procedure  to  be  evaluated.  
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SECTION  II 


DEVELOPMENT  OF  THE  ESTIMATOR 


2.1  INTRODUCTION 

The  object  of  this  chapter  is  to  develop  the  mathematical  models 
needed  to  estimate  the  parameters  of  a synchronous  generator.  The  approach 
taken  is  to  model  the  generator,,  use  this  model  to  predict  outputs  based  on 
the  current  parameter  estimates,  and  then  to  adjust  the  parameter  estimates 
systematically  to  minimize  the  weighted  square  of  the  difference  between 
the  measured  outputs  of  the  generator  and  the  model  outputs.  This  is  the 
essence  of  a weighted  least-squares  estimation  algorithm.  Furthermore, 
incorporating  the  statistics  of  the  noise  which  inevitably  corrupts  the 
measurements  into  the  estimation  algorithm  leads  to  the  maximum  a posteriori 
probability  estimator. 

This  chapter  is  divided  into  two  main  sections.  The  first  presents 
the  model  of  the  generator  in  a form  suitable  for  digital  computer  simula- 
tion. These  equations  are  cast  into  a state-space  notation  for  convenience. 
The  second  section  develops  weighted  least-squares  estimation  algorithm  by 
the  method  of  quasilinearization.  The  statistics  of  the  noise  are  used  to 
derive  an  optimal  estimator,  which  is  the  maximum  a posteriori  probability 
estimator.  This  formulation  is  a special  case  of  weighted  least-squares 
estimation  with  the  weighting  matrices  determined  from  noise  statistics 
and  a priori  information  about  parameter  error  statistics. 

2.2  GENERATOR  MODEL 

A typical  three-phase,  alternating-current  generator  consists  of 
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three  armature  windings  placed  symmetrically  on  the  stator  surrounding  a 
rotor  that  is  driven  externally.  The  rotor  consists  of  an  iron  core  with 
a field  winding  excited  by  direct  current.  Additional  rotor  circuits 
called  damper  windings,  consisting  of  short-circuited  bars,  are  often  im- 
bedded in  the  rotor  surface.  Such  a machine,  considered  here,  is  drawn 
schematically  in  Figure  la. 

Due  to  the  effect  of  the  iron  rotor,  the  armature  inductances  have 
components  that  vary  with  the  rotor  angle.  By  assuming  symmetry  of  the 
rotor  about  the  pole  axis,  or  direct  axis,  and  about  the  interpole  axis 
(or  quadrature  axis)  and  by  assuming  sinusoidally  distributed  armature 
windings  along  the  air  gap,  the  fundamental  component  of  the  air  gap  flux 
linking  the  armature  is  proportional  to  A + Bcos20  [1],  As  a result,  the 
armature  self  and  mutual  inductances  are 


La  = La0  + Lalcos20 


(3) 


Jab 


[L 


abO 


LajCos2 (9 


(4) 


By  projecting  the  armature  circuit  quantities  onto  the  d-q  coordinate 
frame,  which  is  fixed  in  the  rotor.  Park  [1]  obtained  a set  of  differen- 
tial equations  with  constant  coefficients.  Fictitious  windings  along  the 
d and  q axes  have  constant  inductances  since  the  paths  of  flux  have  con- 
stant permeance.  The  self  inductances  of  the  d-q  model  are 

S • La0  + L.b0  + 1 < 


7a 
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Defining 


L1  A La0  + LabO 


(7) 


and 


L 
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al 
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then  Ljj  and  are  computed  from 


h>  " L1  + l2 


LQ  = L1  " L2 


(9) 

(10) 


The  net  result  of  the  preceding  is  the  circuit  model  illustrated  in 
Figure  lb.  It  is  described  by  a set  of  time-invariant,  linear  differential 
equations  and  a set  of  time-varying  measurement  equations.  These  are  given 
by: 


= -(RL  1 + m(t)  + v(t) 
z(t)  = T(t)L-S(t) 
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(12) 
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and  0 - u)t  + The  measurement  equations  correspond  to  a transfor- 

rn.it  ion  ot  coordinates  from  the  D-Q  reference  t lame  to  the  three-phase 
armature  reference  frame,  denoted  hv  the  subscripts  A,  B,  and  C.  The 
rotor  circuit  quantities,  denoted  In’  subscript  K for  the  field  and  by 
tin*  subscripts  Kit  and  KQ  for  the  damper  windings,  are  invariant  under 
this  transformation. 

These  equations  are  not  tin'  same  as  those  derived  by  i’ark  [1] 
hut  have  been  modified  as  suggested  by  Lewis  l 111. 

The  virtue  of  this  particular  transformation  is  that 
the  mutual  reactances  are  always  reciprocal.  This  is  not  true  of 
Park's  original  equations. 

These  equations  are  valid  for  balanced  three-phase  operation 
of  the  machine.  If  an  unbalanced  load  is  connected,  a zero-sequence 
circuit  must  be  added  to  the  0 and  Q circuits  to  represent  the  armature 
circuits.  The  zero-sequence  equations  are  given  below. 
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The  zero-sequence  variables  are  uncoupled  from  the  rest  of  the  equations, 
and 


dll> 
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dt 
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Since  the  analysis  and  experiments  will  be  carried  out  under  balanced 
conditions,  the  zero-sequence  equation  Is  not  considered  further. 

Generator  parameters  are  usually  expressed  In  per  unit,  that  is 
In  dimensionless  ratios  to  base  values.  The  base  values  are  ordinarily 
chosen  to  be  the  rated  values  to  facilitate  comparison  of  machines  of 
different  sizes  and  ratings.  Although  Equations  (ll)  and  (12) 
are  valid  in  any  consistent  set  of  units,  such  as  the  MRS  system,  a per 
unit  system  (s  used  subsequently  to  simplify  the  computations  and  to  be 
consistent  with  accepted  practices. 

Since  i ■ L and  p » , the  direct  axis  part  of  Equation 

(11)  can  be  rewritten  in  operational  form 
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Dividing  each  row  by  the  corresponding  base  voltage  and  dividing  and 
multiplying  each  column  by  the  corresponding  base  current  yields  the 
equations  in  per  unit. 
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Base  current  and  voltage  for  the  armature  circuits  are  the  rated  values. 
Base  field  current  is  chosen  to  be  the  value  of  current  which  causes 
rated  open-circuit  armature  voltage.  Two  constraints  are  placed  on  the 
three  remaining  arbitrary  rotor  base  quantities, 

''KDB^B  = VFB1FB  = VBiB  ' (19) 

These  relations  imply  that  the  base  power  is  the  same  on  all  circuits 
and  that  the  per  unit  inductance  matrix  is  symmetric.  One  more 
constraint  will  be  placed  on  the  base  quantities  of  the  damper  winding 
circuit.  Choosing 
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in  per  unit.  By  similar  reasoning,  choice  of  v i = v„i_  results 

NyD  K^d  d d 

in  a reciprocal  quadrature  axis  inductance  matrix  and  choice 
LQKQ 

1KqB  **  iB  results  in  Lqkq  = L^q.  ln  Per  unit.  As  a result. 

Equations  (ID  and  (12)  are  valid  with  the  following  per-unit  quantities: 
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Equation  (17)  has  two  more  parameters  than  the  per  unit  equations 
since  the  damper  winding  currents  are  scaled  by  a factor  containing  a 
ratio  of  damper  winding  inductances.  This  choice  of  base  currents 
results  in  a reduction  in  the  number  of  inductance  parameters  and  the 
addition  of  an  equal  number  of  parameters  in  the  base  rotor  quantities. 
The  object  of  this  model  is  to  represent  the  generator  by  a circuit 
equivalent  at  the  terminals.  Since  the  damper  windings  are  inaccessible 
short-circuited  turns  embedded  in  the  rotor,  the  current  in  them  does 
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not  need  to  be  determined  absolutely  to  model  the  generator  at  the  termin- 
als. The  net  result  of  this  particular  per  unit  representation  is  that  two 
parameters  are  arbitrary  (any  two  parameters  may  be  arbitrarily  specif  led)  but 
the  resulting  circuit  is  still  equivalent  at  the  armature  and  f ield  terminals . 

Since  the  experiment  was  conducted  with  the  generator  under  a 
balanced  resistive  load,  as  described  in  Section  3,  this  special  case  is 
now  considered.  The  field  is  excited  by  a regulated  voltage  source  pro- 
viding a constant  v^.  The  components  of  the  armature  voltage  are 
proportional  to  the  corresponding  currents. 


vd ' -\‘d 

Vq  ' AS  <22) 

Therefore,  Equation  (11)  still  holds  if  R^  is  replaced  by  R^ + , and  if 
v=  (0  vf  0 0 0)T. 

For  convenience  in  the  following  derivation,  the  model  Equations 
(11)  and  (12)  are  written  in  more  general  notation. 


(23) 

(24) 

x • x(xQ,y,v,t)  = state  vector  (flux  linkages)  = 


= f(x,y,v)  = F(y)x  + v;  x(0)  = xr 


z(t)  = h(x,y,0o,t) 


where 
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V - 


input  vector  (field  voltage) 


t - time 

0Q  " eo(xQ)  » initial  rotor  angle 

x " x (y)  * initial  state  vector 

o o 

z « measurement  vector  (terminal  currents) 

The  matrix  F » -(RL  ^ + Q)  and  the  parameter  vector  is 

y = (L1  Ldf  Lkd  Lf  Lfkd  L2  Lkq  Rd  Rf  \d  Rkq)  <25) 


where 


Ld  * L1  + L2 

Lq  = L1  - 4 (26) 

These  equations  assume  the  generator  is  in  steady  state  at  time 

t=0.  Therefore,  the  initial  state  vector  x and  the  initial  rotor 

o 

angle  0q  are  computed  from  the  parameters  by  the  relations 


x(0)  = XQ(y)  = -F (y)  *v 


(27) 


and 
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For  t>0,  a transient  is  induced  by  a step  change  in  load  resistance 
and  the  state  vector  is  found  by  numerical  solution  of  Equation  (23). 


t 
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Equations  (23)  and  (24)  describe  u continuous-time, 
deterministic  state-space  model  of  the  generator.  Since  a digital 
computer  estimation  algorithm  Is  developed  In  the  sequel,  consider 


the  measurement  Equation  (24)  to  be  a function  of  t^  for 
k«l ,2  , . . . ,k(- , a set  of  discrete  time  points.  Finally,  since  inaccura- 
cies and  noise  inevitably  corrupt  the  measurement  process,  a term 
representing  an  additive  measurement  noise  sequence  is  included. 

The  result  is  expressed  as 


z(t,)“h(x,y,0  , t.  ) + w.  , k - 1,2,  . . . ,kr 
k o k k I 


where  t^  is  the  k discrete  time  point,  and  Is  the  k{^  discrete 
noise  sample. 


2.  ) We  ij^hted  Least  -Squares  Est  lmat  Ion 

Consider  a parameter  estimator,  an  algorithm  not  yet  specified. 


which  recursively  generates  an  estimate  of  the  parameter  vector.  It 
the  current  parameter  estimate  is  denoted  by  y,  then  the  output  from  a 
model  of  the  generator  is  h(x,y,0  ,t,),  where  x is  the  state  vector 

O K 

estimate  computed  from  numerical  solution  of  the  model  Equation  (21) 
using  y.  If  the  measured  output  of  the  actual  generator  is  z(t^), 
then  the  weighted  least -squares  error  criterion  is  given  by 


J - \ )’  { l as  (t.  )-h(x,y,0  , t.  ) ]lQ[z(t  )-h(x,y,0  ,t.)  ] ) . 


Here  Q is  a non-negative  definite  weighting  matrix.  The  interpretation 
of  this  error  criterion  is  that  minimizing  .1  also  minimizes  the  weighted 


squares  ot'  the  error  between  the  observed  output  and  the  modeled  output. 

In  the  estimation  algorithm,  the  method  of  quasi  1 inear 1 zut Ion 
112,13]  will  be  applied.  l.et  the  current  parameter  estimate  y be 
updated  to  produce  the  new  estimate 


y • y + Av 
new 


(31) 


Expanding  Equation  (29)  in  a Taylor  series  about  y and  retaining  only 
the  first  two  terms  gives  the  linear  approximation 
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where 


hk  - h(x,y,0o,tk) 


and 
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Repeated  application  of  the  chain  rule  for  partial  derivatives  give 


dh  ^ 3h  , <>x  ^Xo,  3h  *^o  ^Xo 

dy  3y  3x^3y  3x  dy  36  dx  dy 
7 J o ooy 


(33) 


Likewise,  the  right  hand  side  of  Equation  (23),  f(x,y,v)  may  be  expanded 
Differentiating  Equation  (23)  with  respect  to  v and  x yields  the 
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sensitivity  matrix  differential  equations  below, 
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The  solutions  to  Equations  (34)  and  (35)  are  the  sensitivity 
matrices 


3x  , 3x 
r—  and  r — 
3y  3x„ 


(36) 


A necessary  condition  for  minimizing  J is  that  its  gradient 
with  respect  to  y vanish.  Evaluating  the  gradient  at  the  new  parameter 
estimate  y + Ay  results  in 
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Solving  for  Ay  gives 
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(38) 


Setting  the  new  estimate  equal  to  y + Ay  completes  the  recursive 
algorithm  for  computing  the  parameter  estimates. 

In  summary,  the  algorithm  consists  of  solving  the  differential 
Equations  (23),  (34),  and  (35)  numerically,  computing  sensitivity 
matrices  and  solving  the  system  of  1 inear  algebraic  Equat ions  (38). 
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This  recursive  algorithm  is  ideally  suited  for  implementation  on  a 
digital  computer. 

2 . 4 Statistical  Estimation  Theory 

The  least-squares  approach  of  the  previous  section  ignores  the 
statistical  nature  of  the  estimation  problem.  If  information  about 
the  statistics  of  the  noise  is  available,  the  estimation  process  can 
be  improved.  This  section  presents  a derivation  and  discussion  of 
relevant  aspects  of  stochastic  parameter  estimation  theory. 

2.4.1  Bayesian  Approach 

Consider  an  error  defined  to  be 


y - y(z) 


(39) 


1 


where  y(z)  is  some  estimate  based  on  measuring  z.  Let  C(y-y(z))  be 
the  cost  function  describing  the  penalty  for  making  that  error.  The 
Bayesian  risk  [14,15]  is  defined  as  the  conditional  mean  of  the  error 


VK.  - E{C(y  - y)  | z } 


C(y  - y)p(y|z)dy  . 


For  the  case  of  the  uniform  cost  function. 


C(y-y) 


if  | | y-y 


0 otherwise 


(40) 


(41) 


Consider  the  limit  of  the  cost  function  as  t -O , 
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(42) 


N 

c(y  - y)  - n 6(y  - y ) . 

j«i  J J 


As  a result 


BR  = -p(y(z) |z)  . (43) 

Minimizing  the  Bayesian  risk  is  equivalent  to  maximizing  the  posterior 
probability  density. 

2.4.2  Maximum  A Posteriori  Probability  Estimator 

Denote  the  set  of  all  measurement  vectors  in  the  sequence  by 
Z = (z(t^) .zCt^) , . • . ,z(t^)}.  If  the  measurement  noise  is  represented 
as  a vector  stochastic  process  and  the  parameter  vector  as  a random 
vector,  the  posterior  probability  density  is  given  by  Bayes'  rule  as 

p(y|z)  - P<ZI^P(y)  . (44) 

Since  p(Z)  is  not  dependent  on  y,  maximizing  (44)  by  choice  of  y 
is  equivalent  to  maximizing  p(z|y)p(y). 

If  the  measurement  noise  is  normally  distributed  with  zero  mean 
and  covariance  matrix  V , the  conditional  probability  of  the  k*"*' 
measurement  vector  is  given  by  [14] 

P(z(tk)|y]  = [(2x)mdet(Vw)]“i 

• exp{-i[z(tk)  ~ h(x,y,tk)]TVw1[z(tk)-h(x,y,tk)])  } 

(45) 
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If  the  noise  samples  are  uncorrelated  with  each  other,  that  is  from  a 
white  noise  sequence,  the  probability  of  observing  the  entire  sequence 
Z is  simply  the  product  of  kf  terms  like  the  right  side  of  Equation 
(45). 


plZ|y]  = [(2w)mdet(Vw)] 


-kf/2 


exp{-i  l [z(tk)-h(x,y,tR) ]Tv"1[z(tk)-h(x,y,tk)]} 
k=l 


(46) 


If  the  random  parameter  vector  y is  now  assumed  to  be  normally 
distributed  with  mean  My  and  variance  Vy,  then 

p(y)  = [(2ir)m|Vy|]_iexpt-i(y-My)V1(y-My)]  . (47) 


Since  the  exponential  function  is  monotonic  in  its  argument  and 
since  the  functions  multiplying  the  exponentials  in  Equation  (46) 
and  (47)  do  not  contain  y,  maximizing  p(z|y)p(y)  is  equivalent  to 
minimizing  the  error  criterion 
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J*|  [ [z(tk)-h(x,y,tk)]V1[z(tk)-h(x,y,tk)] 


k=l 


+ |(y-My)V1(y-My)  , 


(48) 


resulting  in 
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(49) 


Thus  the  maximum  a posteriori  probability  estimator  can  be 
considered  a special  case  of  a least-squares  estimator,  including  prior 
information  about  the  parameter  vector,  if  the  weighting  matrices  are 
chosen  equal  to  the  inverse  of  the  error  covariance  matrices. 

2.4.3  Lower  Bound  on  the  Error  Covariance 

A lower  bound  on  the  error  covariance  of  the  estimator  derived 
is  easily  obtained  from  a generalization  of  Fisher’s  information 
matrix  [15,16]. 


E[ (y-y) (v— v )T 1 ^ ( I [ 


k=l 


dhk  _i  dhk  -1 
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dy  w dy  y 


(50) 


The  derivation  of  this  bound  is  presented  in  Appendix  B.  This  lower 
bound  is  computationally  Inexpensive  since  the  right  side  of  Equation 
(50)  is  already  computed  In  the  estimation  In  the  algorithm. 

Equation  (49). 

2 . 5 Summary 

The  derivation  of  an  algorithm  suitable  for  estimating  synch- 
ronous generator  parameters  is  approached  from  the  point  of  view  of 
weighted  least-squares  estimation  theory.  First,  the  equations  des- 
cribing the  generator  are  cast  into  a state-variable  form.  Next, 
the  method  of  quasi l Inear izat ion  is  used  to  derive  a least-squares 
estimation  algorithm.  By  considering  a stochastic  formulation,  the 
maximum  a posteriori  probability  estimator  is  shown  to  be  a special 


case  of  the  weighted  least-squares  algorithm  with  correct  choice  of 
the  weighting  matrices. 

While  the  direct  implementation  of  the  weighted  least-squares 
algorithm  results  from  weights  chosen  arbitrarily  or  by  physical 
intuition,  the  optimal  estimator  improves  the  quality  of  the  estimates 
by  using  the  prior  error  statistics  to  choose  the  proper  weights. 
Unfortunately,  the  exact  statistics  of  the  noise  are  not  known.  The 
approach  taken  here  is  to  study  the  effect  of  incorrect  prior 
statistics  using  a computer  simulation.  Then,  when  processing  data 
from  the  actual  generator,  estimates  of  the  relative  magnitudes  of 
error  variances  can  be  made  to  enable  intelligent  choice  of  weights. 

If  these  estimates  of  the  prior  statistics  are  correct,  the  result 
should  approach  the  performance  of  the  optimal  estimator.  If  the 
statistics  are  in  error,  at  least  the  resulting  suboptimal  estimate 
satisfies  the  least-squares  error  criterion.  In  any  case,  the 
algorithm  minimizes  the  weighted  square  of  the  output  error. 


SECTION  III 


DESIGN  OK  THE  EXPERIMENT 


3. 1 Int  roduct Ion 

The  success  of  an  iterative  parameter  est  imation  algorithm,  such  as 
the  one  previously  described,  depends  on  the  quality  of  the  input  data.  In 
other  words,  the  criterion  of  goodness  for  the  experiment  that  produces  this 
data  must  be  the  performance  of  the  estimator.  In  this  research  the  experi- 
ment was  designed  with  the  aid  of  computer  simulat  ions.  This  enabled  various 
experimental  conditions  to  be  tested  with  maximum  flexibility. 

First,  certain  constraints  on  available  equipment  and  on  eomput  tng  and 
processing  t ime  were  recognized.  These  were  primari  !y  economic  l imitations. 

The  experiment  was  designed  within  these  constraints  by  selecting  such  factors 
as  machine  load  and  excltat  ion,  data  sampl  ing  rate  and  overal  l data  record  length. 
To  choose  these  condit  ions  Intel  l igent  Iv,  a computer  simulat  ion  of  a typical 
generator  was  used  to  model  the  experimental  setup.  The  result  ing  test  data 
were  used  to  assess  the  performance  of  the  est  imator.  Therefore,  the  effect  of 
t he  undetermined  exper  i ment  a 1 condi  t ion  was  eas  i 1 y est  ab  l ished  and  t he 
experiment  thereby  designed. 

1.2  Equipment  Constraints 

The  test  generator  was  a three-phase,  four-pole,  alternat  ing  current . 
synchronous  machine  rated  at  3 KVA  and  230  volts  at  40  hertz.  The  drive  motor  was 
a IS  horsepower  direct  current  machine  rated  at  1200  revolut  ions  per  minute.  The 
do  supply,  froma  large  motor-generator  set  , had  only  open-loop  control . As  a 
result , the  speed  of  the  drive  motor  was  manual  Iv  adjust  able  with  no  provision 


for  automatic  regulation. 


The  load  bank  consisted  ot  three  power  resistors  connected  in 


wye  with  an  additional  resistor  suddenly  switched  in  parallel  with 

each  leu  to  induce  a transient.  This  arrangement  is  Illustrated  in 

Figure  2.  The  closing  of  the  three-pole  switch,  labeled  S,  causes 

R1R2 

the  load  resistance  to  suddenly  decrease  from  R.+R,  to  - — + R,. 

I 3 K j tK  , J 

Resistors  R^  are  current  shunts  used  to  obtain  voltage  signals  pro- 
portional to  the  load  current  in  each  leg.  These  signals,  along  with 
a similar  one  proportional  to  field  current,  were  recorded  and  processed 
as  described  in  the  next  section. 

Excitation  of  the  field  was  supplied  by  a regulated  electronic 
dc  power  supply  rated  at  0.5  amp  and  500  volts.  This  was  adequate  to 
supply  up  to  the  rated  field  current  of  .525  amps  at  the  rated  field 
voltage  ot  125  volts. 

The  main  constraint  posed  by  the  available  equipment  was  to 
limit  armature  current  to  values  within  the  generator  ratings  to  pre- 
vent damage  to  the  windings.  The  second  constraint  was  a limit  on 
the  maximum  change  in  armature  current  during  the  test.  To  insure 
that  these  constraints  were  met,  the  field  voltage  was  reduced  below 
its  rated  value,  reducing  the  magnitude  of  the  load  currents. 

3.3  Simulation  I'sing  Nominal  Parameters 

Within  the  constraints  already  posed,  several  experimental 
conditions  remain  to  be  selected.  First,  the  magnitude  of  the  step 
load  change  must  be  selected  large  enough  to  enable  all  the  parameters 
to  be  estimated.  Some  parameters,  notably  the  damper  circuit  rectances 
and  resistances,  affect  the  terminal  currents  only  during  transients. 
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Thus  a significant  transient  must  be  introduced  by  the  sudden  load 
changes.  Second,  the  data  sampling  rate  must  be  chosen  fast  enough  to 
accurately  represent  the  terminal  current  waveforms.  The  well-known 
sampling  theorem  states  that  the  minimum  rate  is  twice  the  highest 
frequency  component  of  the  waveform.  Experience  shows  that  the  practi- 
cal minimum  is  somewhat  faster  than  the  theoretical  limit.  A balancing 
constraint  is  the  limitation  of  the  total  amount  of  data  storage  space 
available.  Finally,  the  length  of  the  data  record  and  the  number  of 
load  switchings  must  be  determined. 

To  determine  these  conditions- so  that  the  experiment  would  be 
successful,  the  experimental  setup  was  modeled  with  a computer  simula- 
tion. The  machine  was  modeled  with  Equations  (23)  through  (29) 
solved  numerically  using  a modified  Runge-Kutta  method  [17].  The  modi- 
fication, due  to  Merson  [18],  allowed  an  estimate  of  the  roundoff  error 
to  be  computed  to  insure  the  step  size  was  small  enough.  A multiplicative 
type  pseudo- random  noise  generator  [19]  simulated  additive  measurement 
noise.  The  simulation,  implemented  on  a small  digital  computer  with  a 
magnetic  disk  operating  system,  provided  test  outputs  to  allow  assess- 
ment of  the  effect  of  the  experimental  conditions  on  the  parameter 
estimator.  This  information  allowed  the  experiment  to  be  designed. 

To  implement  the  numerical  solution  of  the  machine  model,  a set 
of  typical  parameters  was  determined  from  nameplate  data,  a few  simple 
tests,  and  a list  of  typical  machine  constants  [3].  The  parameters 
Xj,  Xjf*  Rj  and  were  measured  by  simple  ac  and  do  steady-state 
tests.  The  remaining  parameters  were  roughly  estimated  from  the 
typical  parameter  list.  The  results  of  this  nominal  parameter 

30 


1 


computation  are  given  in  Table  1.  The  results  of  the  steady-state  tests  and 
the  details  of  the  calculation  are  presented  in  Appendix  C.  It  should  be 
emphasized  that  the  nominal-parameter  model  of  the  machine  is  not  an  accu- 
rate representation  of  this  generator  but  is  similar  to  a typical  generator. 

The  first  run  of  the  simulation  modeled  the  generator  during  a sudden 
three-phase  short  c ircuit  on  the  armature  terminal s . The  main  purpose  of  this 
run  was  to  test  the  operation  of  the  estimator.  No  measurement  noise  was  added 
and  the  measurements  were  weighted  equal  ly.  The  step  size  was  2 msec,  and  the 
field  voltage  was  the  rated  value.  The  parameter  estimates  converged  rapidly 
as  shown  in  Figure  3a  and  3b,  which  are  plots  of  the  computed  relative  error 
versus  iteration.  The  relative  error  is  defined  as 

f - ^ • (51) 

where  y Is  the  parameter  estimate,  and  y is  the  actual  parameter  value,  1 .e. 
the  parameter  value  of  the  model  simulation.  Figure  3c  shows  the  instantaneous 
phase  A current  and  the  field  current  using  parameter  values  which  are  est  i mated. 
The  outputs  based  on  the  parameter  est  i mates  match  the  data  within  graphical 
error. 

The  second  test  is  similar  to  the  first;  however,  a transient 
was  Induced  by  switching  load  resistance  rather  than  a sudden  short  cir- 
cuit. The  switched  load  test  is  a more  realistic  simulation  of  the 
actual  experimental  setup.  With  the  machine  initially  in  steady  state, 
the  load  resistance  was  suddenly  switched  from  1.0  to  0.23  per  unit. 

The  toad  was  switched  back  to  one  per  unit,  400  milliseconds  later. 
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TABLE  1 


NOMINAL  PARAMETERS  IN  PER  UNIT 


Parameter 

Nominal  Value 

Number 

Parameter 

(pier  unit) 

1 

L1 

2.69xl0~2 

2 

Ldf 

4.87xl0~3 

3 

k)cd 

2.2xl0-3 

4 

Lf 

10. 8xl0_3 

5 

Lfkd 

4.01xl0“3 

6 

L2 

. 565xl0~3 

7 

hcq 

1.17xl0"3 

8 

Rd 

1.42xl0-2 

9 

Rf 

2.31xl0“2 

10 

\d 

7 . 59xl0~2 

11 

% 

7.59xl0“2 
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Convergence  of  the  estimator  was  somewhat  slower  than  in  the 
previous  test,  as  shown  in  Figure  4a  and  4b.  Figure  4c 
shows  the  terminal  currents  computed  with  the  parameter  estimates  after 
ten  Iterations  to  be  quite  close  to  the  correct  values.  The  errors 
remaining  in  the  estimates,  then,  have  a small  effect  on  the  terminal 
currents. 

Tests  run  with  switched  loads  of  differing  magnitudes  show  that 
smaller  resistances,  approaching  the  short-circuit  test  in  the  limit, 
produce  results,  closer  to  those  obtained  in  test  one.  On  the  other 
hand,  resistances  approaching  1.0  lead  to  smaller  transients  and 
poorer  results. 

Tests  in  which  one  transient  is  Induced  by  load  switching  show 
poorer  results  than  the  previous  test  with  two  transients.  Two  load 
switchings,  from  normal  load  to  small  load  to  normal  load,  provide 
redundant  informat  ion.  This  aids  the  operation  of  the  estimator, 
particularly  in  the  presence  of  noise.  More  than  two  load  switchings 
are  not  considered,  because  more  than  two  operations  of  the 
three-pole  switch  are  not  practical.  The  experiment  then  was 
constrained  to  two  load  resistance  switchings. 

Further  tests  with  different  switching  times  showed  essentially 
the  same  behavior.  The  minimum  time  for  operating  the  switch  prohibited 
very  short  durations,  and  the  limitations  on  overall  data  record  length 
prohibited  very  long  durations.  For  manual  operation  of  the  switch,  a 
duration  of  several  hundred  milliseconds  is  reasonable. 

The  next  test,  illustrated  in  Figure  So  shows  the  effect  of 
ten  percent  errors  in  the  initial  parameter  guesses  and  measurement 
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Figure  4c. 
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noise  with  zero  mean  and  standard  deviation  1 x 10  . Figures  5c  and  5d 

show  instantaneous  test  values  of  a phase  A current  and  field  current, 
respectively.  All  measurements  and  parameters  were  weighted  equally. 

The  effect  of  increasing  measurement  noise  is  to  degrade  the 
performance  of  the  estimator.  At  noise  of  standard  deviation  greater 
than  1.0,  the  estimator  does  not  converge.  With  noise  of  standard 
deviation  of  0.1,  the  estimator  converged  very  slowly,  but  large 
errors  in  some  parameters  show  that  they  tend  to  track  the  noise. 

The  final  test  with  simulated  data,  shown  in  Figure  6a,  was  the 
same  as  the  preceding  test;  however,  the  weighting  matrices  were  chosen 
as  the  inverse  covariance  matrices.  Figures  6b  and  6c  show  instantaneous 
test  values  of  phase  A current  and  field  current,  respectively.  This 
is  an  implementation  of  the  maximum  a posteriori  probability  estimator. 

The  results  show  only  minor  improvements  on  the  preceding  weighted 
least-squares  approach,  shown  in  Figures  5a-5c. 

The  results  of  the  simulated  experiment  show  that  the  estimator 
will  work  satisfactorily  with  a resistive  load  switched  from  1.0  to 
0.25  to  1.0  per  unit,  with  practical  switching  periods,  in  the  presence 

f 

of  moderate  initial  parameter  errors  and  measurement  noise.  This 
experiment  was  consistent  with  the  limitations  imposed  by  available 
equipment.  The  implementation  of  this  experiment  is  discussed  in  the 
next  section. 
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IMPLEMENTATION  OF  THE  EXPERIMENT 

4. 1 Introduction 

This  section  describes  an  experiment  to  measure  the  parameters  of 
a synchronous  generator  in  the  laboratory.  Terminal  current  measurements 
are  recorded,  stored  in  digitized  form,  and  fed  off-line  into  the  parameter 
estimation  algorithm.  The  resulting  parameter  estimates  are  used  to  pre- 
dict new  results,  which  are  compared  to  actual  measurements  to  validate 
the  model. 

The  three  main  thrusts  of  this  part  of  the  research  are  reported 
in  the  remainder  of  this  section.  First  the  data  collection,  processing 
and  storage  methods  are  described  in  detail.  Next,  the  results  of  the 
experiment,  the  parameter  estimates,  are  presented.  Finally,  the  ade- 
quacy of  the  results  is  assessed  by  using  the  parameter  est invites  in  a 
computer  simulation  to  predict  the  generator  response  to  a relatively 
large  change  in  load.  This  prediction  is  compared  to  actual  measurements 
under  the  same  conditions. 

4 . 2 Data  Processing  Methods 

Analog  signals  proportional  to  the  three  armature  currents  and  the 
field  current  are  produced  by  current  shunts  in  the  four  generator  terminal 
circuits.  These  signals  are  recorded  on  separate  channels  of  an  FM  instru- 
mentation tape  recorder.  A typical  channel  is  shown  schemat leal  1 y in  Figure 
7. 

Since  the  parameter  estimator  is  implemented  on  a digital  computer, 
the  four  channels  of  analog  data  are  digitized  off-line  and  stored  in 
serial  form  on  a magnetic  disk  cartridge.  Figure  8 is  a schematic  diagram 
of  the  analog-to-d  igital  (A/D)  conversion  and  the  multiplexing  of  the 
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Figure  7.  Schematic  Diagram  of  Instrumentation  for  Channel  i 


four  digital  data  channels  into  one  serial  record  on  the  disk.  The  A/D 
converters  provide  sampled-data  output  in  12  bit  two's  complement  binary 
numbers.  The  A/D  driver  software  writes  these  numbers  to  the  disk 
serially.  The  stream  of  data  to  disk  contains  sample  one  of  channel  one, 
then  sample  one  of  channel  two  and  so  forth.  Sample  two  of  channel  one 


follows  sample  one  from  the  last  channel.  This  process  of  multiplexing 
data  to  the  disk  is  illustrated  for  two  channels  of  triangular  waveform 
data  in  Figure  9. 

Thus  the  data  stored  on  the  disk  forms  a time  series: 


/x  . ^ AT.  . . 2AT.  . , 3 AT 

^1  * z2 (Ck  4^'  z3 (Ck  4 ^ * ^4  ^k  ^ 4 ' 


zl(tk+l)’  z2(tk+l  + 4} z4(tf)  * 


Here  AT  is  the  overall  time  step  and  = tk  + nata  this 

serial  multiplexed  format  is  used  directly  in  a weighted  least-squares 
algorithm  by  updating  the  state  vector  and  the  sensitivity  matrices  at 
all  the  time  points 


, AT  ^AT  JlAT 

V k + 4 ’ k 4 ’ tk  4 ’ k+1 f • 


The  error  between  observed  and  computed  outputs  is  defined  as 
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Figure  9.  Reduction  of  Two  Channels  of  Digitized  Data  to  One 
Serial  Record,  (a)  Two  Channels  of  Data  Sampled 
Alternately,  (b)  Multiplexed  Serial  Data. 
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(52) 

I 

This  new  definition  is  used  to  implement  the  estimator  directly  from 
the  data  on  disk. 

4. 3 Instrumentation 

The  instrumentation  used  in  this  research  includes  an  ac  ammeter, 
an  ac  voltmeter,  and  a digital  multimeter  for  dc  measurements.  These 

■ 

meters  are  listed  with  model  numbers  and  manufacturers  in  Table  2.  In 
addition  to  these  instruments  for  steady-state  measurements,  a Honeywell 
model  5600  instrumentation  FM  tape  recorder  was  used  to  record  transient 
data.  Analog-to-digital  conversion  and  data  processing,  described 
previously,  were  implemented  on  a Data  General  NOVA  small  computer. 

4.4  Result  of  the  Experiment  ! 

I 

After  operating  the  generator  at  normal  load  for  fifteen  minutes 
to  minimize  variations  due  to  temperature  changes,  the  load  resistance 
was  measured  as  0.98  per  unit.  The  additional  load  resistor  bank  was 
switched  in  parallel  to  this  load,  and  the  resistance  of  the  combination 
measured  as  0.23  per  unit.  Then  the  load  was  returned  to  normal,  and 
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TABLE  2 

LIST  OF  INSTRUMENTS  USED 


DESCRIPTION 

NUMBER 

MANUFACTURER 

AC 

ammeter 

AA-10 

General 

Electric 

AC 

voltmeter 

AV-12 

General 

Electric 

Digital  multimeter 

3476B 

Hewlett- 

-Packard 

5.1 


the  rotor  speed  measured  as  252  radians  per  second.  The  first  test  was 


performed  under  these  conditions.  The  resulting  data,  recorded  and 
processed  as  previously  described,  was  sampled  at  a rate  of  500  Hertz 
per  channel  or  an  overall  rate  of  2KHz  for  all  four  channels. 

To  determine  a priori  information  about  the  measurement  noise, 
the  recorder  was  run  with  the  inputs  shorted  to  ground.  After  digiti- 
zing, this  data  indicated  approximately  equal  noise  variances  in  all 
four  channels.  As  a result  the  weights  for  the  measurements  were 
selected  as  unity.  Due  to  lack  of  good  a priori  information  on  the 
initial  parameter  error  variances,  these  weights  were  adjusted 
empirically  to  obtain  good  results  from  the  estimator.  By  weighting 
the  parameters  heavily  at  first,  estimates  were  prevented  from  large 
excursions  on  the  first  few  iterations,  enhancing  the  stability  of  the 
estimator.  The  weights  chosen  are  summarized  in  Table  3,  while  the 
parameter  estimates  are  plotted  versus  iteration  number  in  Figure 
10a. 

Figure  Ha  shows  a comparison  of  the  phase  A current  and  the 
field  current  predicted  from  the  parameter  estimates  to  the  corresponding 

i 

data  from  the  experiment.  First,  in  Figure  11a,  the  results  of  ' 

the  initial  parameter  estimates  are  compared  to  the  experimental  data.  ' j 

j 

Then  Figure  lib  shows  the  results  of  the  final  parameter  estimates 

» 

compared  to  the  data.  The  final  results  follow  the  data  quite  well. 

This  indicates  that  the  model  fits  the  data  well  at  these  particular 
operating  conditions. 

To  test  these  results  under  another  condition  and  to  assess  the 
validity  of  the  model,  a second  test  was  run  with  load  resistance 


54 


TABLE  3 

PARAMETER  ERROR  WEIGHTS 


ITERATION 

PARAMETERS 

WEIGHT 

1-13 

Y -Y 

1 7 

1x10 

P 

1-3 

VYn 

1x10 

4-13 

VY9 

1x10 

4-13 

Y10'Y11 

1x10 

ti 


ss 


8 


9 1 


1 


I (Iteration  Number) 

Figure  10a.  Parameter  Estimates  Versus  Iteration  Number  from 
Experimental  Data  (Test  One),  Parameters  1,  2,  3 
5,  6,  and  7 
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I (Iteration  Number) 

Figure  10b.  Parameter  Estimates  Versus  Iteration  Number  from 
Experimental  Data  (Test  One),  Parameters  4,  8,  9, 
10  and  11 


Figure  lib.  IMwjso  A and  Field  Currents  for  Tost  One, 
Iteration  13 


switched  from  0.99  per  unit  to  0.175  per  unit.  This  test  was  first 
simulated  on  the  digital  computer  using  the  parameter  estimates  from 
the  first  test,  then  implemented  in  the  laboratory.  Figure  12 
shows  the  comparison  of  the  simulated  results  to  the  experimental 
results.  The  reasonable  agreement  of  these  results  indicates  the 
success  of  the  model  of  the  synchronous  generator  at  similar  operating 
conditions.  This  does  not,  however,  guarantee  validity  of  this  model 
over  much  larger  changes  in  operating  conditions.  In  fact,  a synch- 
ronous generator  exhibits  nonlinearities;  therefore,  this  linearized 
representation  is  probably  somewhat  inaccurate  for  very  large 
perturbations.  The  model  is  valid,  though,  about  the  conditions  of 
the  tests. 

II 
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SECTION  V 


CONCLUSIONS 


5. 1 Summary  of  Results 

This  work  is  a new  approach  to  an  old  problem:  system  identifi- 
cation theory  applied  to  the  experimental  determination  of  synchronous 
generator  parameters.  The  power  of  this  statistical  identification 
technique  allows  the  parameters  of  Park's  equations,  cast  in  a state- 
variable  formulation,  to  be  determined  using  current  measurements  of 
the  three  armature  phase  currents  and  the  field  current. 

Considering  the  identification  problem  as  a multiport  boundary- 
value  problem  and  applying  the  method  of  quasilinearization  leads  to  a 
weighted  least-squares  parameter  estimator.  Including  additive  noise 
corrupting  the  measurement  process,  a statistical  analysis  proves  the 
maximum  a posteriori  probability  estimator  is  achieved  by  selecting  the 
weighting  matrices  as  the  inverse  noise  and  parameter  error  covariance 
matrices.  This  estimator  is  the  optimal  estimator  in  the  sense  of 
minimizing  the  Bayesian  risk.  If  the  error  covariances  are  incorrect, 
the  estimator  is  no  longer  optimal  but  is  still  a good  implementation 
of  a least-squares  estimator. 

Before  the  experiment  was  implemented,  several  experimental 
conditions,  such  as  the  magnitude  and  duration  of  the  transient  and  the 
data  sampling  rate  and  duration,  needed  to  be  selected.  Also,  a set  of 
constraints,  mainly  economic  in  nature,  were  recognized.  Within  these 
constraints  the  experiment  was  designed  with  the  aid  of  a digital  com- 
puter simulation.  That  is,  a numerical  solution  of  the  mathematical 
model  of  a typical  synchronous  generator  under  a switched  resistive 
load,  was  used  to  test  various  experimental  conditions.  The  results 
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show  that  a sufficient  amount  of  data  is  produced  by  a load  resistance 
switched  from  full  load  value  to  0.25  per  unit  then  switched  back  to 
full  load  value.  The  duration  of  the  transient  was  consistent  with 
manual  operation  of  a three-pole  switch.  The  resulting  experiment 
provided  sufficient  data  at  a sampling  rate  which  could  be  handled  by 
the  data  recording  and  processing  facilities  available,  while  being 
implemented  on  available  laboratory  equipment. 

Implementing  this  experiment  in  the  laboratory,  recording  the 
terminal  current  data,  and  digitizing  this  data  for  input  to  the  off- 
line parameter  estimator  led  to  estimates  of  the  parameters  of  the 
mathematical  model  of  the  generator.  These  parameter  estimates  provide 
a linearized  model  of  the  generator  which  is  valid  over  a range  of 
operating  conditions  near  the  experimental  conditions.  The  results  of 
predicting  the  generator  response  to  a larger  step  change  in  load 
resistance  compare  favorably  to  the  actual  measured  response  under 
those  conditions.  This  favorable  comparison  validates  the  results, 
showing  that  the  model  can  indeed  predict  the  generator  behavior  under 
similar  load  conditions. 

5.2  Significance  of  Results 

In  fact,  the  synchronous  generator  is  a nonlinear  device,  exhibiting 
saturation  and  other  deviations  from  the  assumed  linear  model.  By 
assuming  this  linear  model  structure  and  using  an  iterative  estimation 
algorithm  to  calculate  parameters,  the  results  are  effectively  lineari- 
zed about  the  test  conditions.  That  is,  the  estimator  fits  the  best 
linear  model  to  the  data.  This,  of  course,  presents  a limitation  to 
modeling  drastically  different  operating  conditions.  For  relatively 


small  perturbations  in  conditions,  however,  a simple  valid  solution 
has  been  achieved.  It  should  be  emphasized  that  the  approach  used  here 
is  an  improvement  over  methods  which  assume  linearity  of  the  model  over 
large  changes  in  operating  condition.  For  example,  estimates  of  X^,  the 
synchronous  reactance,  can  be  obtained  from  combining  results  of  open- 
circuit,  unsaturated,  armature  voltage  and  of  short-circuit  armature 
current.  However,  such  tests  assume  that  results  of  two  tests  at 
drastically  different  points  can  be  combined  to  model  the  generator 
under  still  different  conditions,  such  as  at  unity  power  factor  and 
full  load.  This  implies  that  the  generator  model  is  linear  to  large 
perturbations.  Thus,  the  present  work  should  be  viewed  as  a step 
toward  proper  modeling  of  a nonlinear  device  by  a linear  model  valid 
over  a small  operating  region. 

The  assumption  of  a linear  model  was  made  primarily  for  simplicity. 
The  estimation  of  parameters  of  a dynamical  system,  such  as  the  generator, 
can  be  viewed  as  an  inherently  nonlinear  problem,  even  when  the  model 
itself  is  linear.  Augment  the  state  vector  x with  the  parameter 
vector  y. 


x = [--] 
y 


(53) 


Now,  in  the  linear  model  used,  the  parameters  are  multiplied  by  the 
states  in  the  right  side  of  Equation  (23).  As  a result,  this  can 
be  rewritten  as  a nonlinear  dynamical  system: 


dx 


dt 


= f(x) 


(54) 
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The  problem  can  now  be  formulated  as  a nonlinear  state  estimation 
problem.  As  a result  of  this  view  of  the  problem,  any  nonlinearity 
that  can  be  formulated  can  be  included  by  this  approach. 

Therefore,  a general  method  for  determination  of  synchronous 
generator  parameters  from  the  results  of  test  data  has  been  presented. 
The  particular  model  used  here  is  Park's  linear  circuit  model  of  an 
ideal  generator.  A solution,  using  available  equipment,  shows 
favorable  results. 

5 . 3 Recommendations  for  Future  Research 

Inclusion  of  nonlinearities,  especially  saturation,  in  the 
generator  model  deserves  additional  research.  The  preceding  section 
shows  that  nonlinear  differential  equations  are  easily  handled  by  this 
algorithm.  Consequently,  future  research  should  concentrate  on  finding 
a suitable  parameterization  of  the  nonlinearities.  Once  this  is  done, 
the  extension  of  the  present  algorithm  to  include  estimating  the 
parameters  that  describe  the  nonlinearities  is  straightforward. 

In  addition,  the  development  of  algorithms  for  parameter  estimation  for 
permanent  magnet  synchronous  generators  and  brushless  machines  where,  in 
each  case,  the  field  current  is  not  available  as  a measurable  quantity 
are  required. 
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SF.CT  I ON  VI 


ESTIMATION  OF  IEEE  STANDARD  PARAMETERS 


6. 1 Introduction 

For  the  purpose  of  digital  simulation  and  application  of  modern 
systems  theory,  it  is  most  convenient  to  represent  a synchronous  machine 
in  state-space  form.  The  electrical  behavior  of  the  machine  is  then 
determined  by  coefficient  matrices  whole  elements  can  be  expressecf^ in 


terms  of  standard  machine  parameters  as  defined  in  the  IEEE  code. 


21,22 


23,24 

An  output  sensitivity  analysis  study,  ’ of  an  alternator-load 

system  showed  that  its  electrical  output  variables  are  very  sensitive 

to  variations  of  certain  standard  machine  parameters.  Consequently,  the 

set  of  parameter  data  with  which  the  model  has  to  be  provided,  should  be 

sufficiently  accurate.  Procedures  to  determine  the  standard  parameters 

of  a synchronous  machine  have  been  establ ished^ where  several  different 

tests  have  to  be  performed.  With  increasirgly  detailed  model  structures 
25-27 

of  the  machine,  several  methods  were  introduced  to  find  additional 
28-  30 

parameter  values.  In  general,  however,  the  need  for  a method  to 

obtain  more  accurate  parameter  values  for  a given  model  structure  is 
24  25 

recognized.  * One  of  the  primary  tasks  of  an  IEEE  Working  Croup  is  to 
suggest  standard  methods  to  determine  the  parameter  values  for  synchro- 
nous machine  models. ^ 

A newly  developed  test  method  to  accurately  determine  the  coeffi- 
cients of  an  alternator  transfer  function  is  a frequency-response  tech- 
31  32 

nique,  where  judgement  is  required  in  locating  the  breakpoints  of 

the  plot  to  get  a good  fit  of  the  data.  Its  practical  implementation  need 
a variable  frequency  source  capable  of  supplying  relatively  large  currents. 
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In  this  section,  a weighted-least-squares  and  maximum  likelihood 


estimation  technique  is  presented  to  determine  the  complete  set  of 

standard  parameters  of  a salient-pole  synchronous  machine  with  damper 

winding  from  the  observed  output  data  of  a sudden  short-circuit.  An 

advantage  of  identifying  standard  parameters  is  that  the  parameter 

values  obtained  can  be  used  for  any  model  formulation  of  the  same  order 

employing  standard  parameters,  which  are  best  known  to  power  engineers. 

In  this  connection,  it  is  noted  that  the  values  of  the  coefficient- 

matrix  elements  of  the  corresponding  state-space  model  are  uniquely 

determined  from  a given  set  of  values  of  the  standard  parameters  whereas 

24 

this  is  generally  not  true  for  the  reversed  case.  The  chosen  test 
procedure  has  the  feature  that  regular  test  equipment  can  be  employed 
in  the  practical  implementation  and,  above  all,  fast  convergence  of 
the  estimation  scheme  is  achieved  while  observation  time  and  estimation 
error  are  kept  to  a minimum. 


. 


68 


6.2 


Parameter  Estimation 


6.2.1  Weighted-Least-Squares  (WLS)  Method 

Let  a real  physical  system  be  described  by 

x = A(a)x  + B(a)u;  x(tQ)  = xo  , (55) 

y = C («)x 

where  x and  y are  the  state  and  observable  output  vectors,  respectively, 
and  u represents  the  input  vector  assumed  to  be  a given  time  function. 
A(a),  B(a)  and  C(a)  are  coefficient  matrices  dependent  on  a parameter 
vector  a,  whose  value  is  to  be  identified.  Note  that  the  unknown  para- 
meter a does  not  necessarily  correspond  directly  to  the  elements  of  any 
of  the  coefficient  matrices. 

33 

The  value  of  a is  chosen  to  minimize  a weighted-least-squares 

error 

N 

J = I tyr(i)  - y (i) ]TW(i) lyr(i)  - y(i)]  . (5b) 

i=l 

Herein,  yr(i)  and  y(i),  i=l,2,...,N,  represent  the  sampled  outputs  of 
the  real  system  and  model  reference,  respectively,  both  excited  with 
the  same  input,  whereas  W(i)  is  a positive  semidefinite  symmetric 
weighting  matrix  chosen  on  the  basis  of  engineering  judgment.  Expanding 
the  output  y (i ) in  a Taylor  series  about  its  trajectory  at  a given  value 
aQ  of  parameter  a gives 


y(i)  = y(i,aQ)  + [3y(i)/9a]a  (a  - a ) 

o 


+ higher  order  terms. 


where 


[9y(i)/3ot]a  = C (aQ)  [9x  (D/gci]^  + [9c  (a)/9a)a  x (i  ,ao)  (58) 

o o o 


The  time-dependent  variables  [3x(i)/9a]  are  found  by  solving  and 

o 

sampling  the  partial  derivative  of  (55)  relative  to  a,  evaluated  at 


d (9x/9a ) /dt  = A(a  ) (9x/9a)  + [9A(a)/9a]  x(a  ) 

CL  O OL  OL  O 

o o o 


+ [SBCaJ/SaJ^  u;  (3x/3a)a  = 0 , 

o o ' o 


and  equation  (55)  itself,  evaluated  at  aQ.  Substitution  of  (57) 


in  (56)  yields 


J = [ (y_(i)-y(i,a  ) - [9y(i)/9a]  (a  - a ) TW(i) 

. r o a o 


• (yr(i)  - y (i ,aQ)  - [9y(i)/9a]a  (a  - aQ) } , 


where  the  higher-order  terms  of  the  expansion  are  ignored.  With 
(o  - aQ)  selected  such  that  (60)  is  minimized. 


9J/9 (a  - a ) = 0 
o 


I 


from  which  follows  a set  of  linear  algebraic  equations  for  (a  - a ) : 

o 


{ £ [3y(i)/3a]^  W(i)/3a]aQ} (a  - aQ)  = 


i=l 


l [3y(i)/3a]a  W(i)[yr(i)  - y(i,aQ)] 
i=l  o 


(62) 


where  a denotes  the  a that  minimizes  (60).  From  (55),  (58), 

(59)  and  (62),  an  iterative  solution  for  a can  be  obtained.  Upon 
writing  (a  - uq)  = Aa,  the  recursion  formula  for  computing  successive 
estimates  of  a is 


“k  ~k-l  *k 
a = a + Aa 


(63) 


where  k(=l,2,...)  refers  to  the  iteration  count  and  a is  some  initial 
guess  made  by  engineering  judgment.  Since  the  derivation  assumes  a 
sufficiently  small  Aa  due  to  truncating  the  Taylor  series,  it  is  pro- 
posed to  modify  (63)  in 


*k  'k-1  ^ „„  ,'k 
a = a + G(k)Aa 


(64) 


where  G(k)  is  a gain  matrix,  whose  elements  are  selected  on  the  basis 
of  computational  judgment.  G(k)  controls  the  extent  of  correction  in 
each  iteration  step.  An  excessively  large  value  of  G(k)  may  cause  a 
to  overshoot  and  go  into  oscillation  about  a,  whereas  an  overly  small 

■'K 

value  of  G(k)  causes  very  slow  convergence  of  a to  a.  One  of  the 
many  possible  variations  is  to  take 
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(65) 


G(k)  = diag(gii(k)] 

where  the  scalars  g^(k)  are  chosen  In  some  proper  manner. 

6.2.2  Maximum  Likelihood  (ML)  Method 

It  was  assumed  in  the  previous  section  that  neither  input  nor 
output  of  the  system  is  corrupted  with  noise.  Considering  the  physical 
nature  of  the  case  to  which  the  estimation  technique  is  to  be  applied, 
ignoring  the  input  noise  seems  to  be  quite  acceptable.  However,  the 
output  measurement  is  generally  noise  corrupted.  Therefore,  a more 
realistic  expression  for  the  system  output  is 

y = C (a) x + £ (66) 

Here,  £ represents  a measurement  noise  which  is  assumed  to  be  zero-mean 
gauss ian  white  with 

E{£(t)£T(T) } = R(t)5(t  - T)  , (67) 

where  R(t)  is  positive  definite,  and  E{*}  and  S(t-x)  stand  for  the 
expectation  operator  and  Dirac  delta  function,  respectively. 

The  value  of  a is  chosen  from  the  observed  sampled  output  yN  = 

T T T T 34 

|y  (l),y  (2),...,y  (N) ) such  that  the  likelihood  function  (l(y  :a), 

which  depends  on  a , is  maximum.  A priori,  the  functional  relationship 

between  y and  a is  given  by  the  conditional  probability  density 

function  p(yN|u).  Thus,  using  the  product  rule  for  probabilities,  the 

likelihood  function 
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(68) 


I 


L(yN:a)  = n p{y  (i) |yi_1,a) 
i=l 


T T T T i 

where  y^  = [y  (l),y  (2),...,y  ( i ) ] and  p{y  (i)  | > is  the  conditional 


density  function  of  y(i)  for  given  y^_^  and  a.  For  the  system  described 
by  (55)  and  (66), 


m I d r i \ 1 1 ~i 


p{y  (i) |yi_1a}  = [ (2u)m| R(i) | ] 


exp{-i [y (i) -y (i) ]TR  1 ty(i)-y(i) ] } 


(69) 


where  the  mean  output 


y(i)  = C(a)x(i)  , 


(70) 


and  m is  the  dimension  of  output  vector  y.  Therefore,  it  follows  from 
(68)  and  (69)  that  maximization  of  L(yN:a)  and  minimization  of 


N 


J = l [y(i)  - y (i) ]TR  1 [y ( i)  - y (i) ] 
i=l 


(71) 


both  with  respect  to  a are  equivalent.  Considering  y = C(a)x  and  (70), 
expression  (71)  is  identical  to  (56)  if 


W = R 


-1 


(72) 


Therefore,  the  same  computational  algorithm  as  for  the  WLS  estimation 
can  be  applied  if  the  inverse  covariance  matrix  of  the  measurement  noise 
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! 1 


! 


determined  from  the  a priori  knowledge  is  vised  foi  the  weighting 


matt i x. 


t> . 3 Apj> l_i ca tion  to  Synchronous  Machine s 
6.3.1  Machine  Model  at  Sudden  Short-Circuit 

A salient-pole  synchronous  machine  with  damper  winding  is 
considered.  In  the  following,  the  machine  equations  are  written  in 


state- space  form  as  required  by  (SS)  with  f lux-1  inkages  selected 


as  the  state  variables 


20 , 2 3 


*d  * ‘W  V*V  (cdkd/Td)  *kd+  (rdf/Td)  ♦f  *vd 


<J’  - -un|)>(c  /T  )i|'  + (c  /T  ) 4'.  tv 

q q q q qkq  q kq  g 


1',  - (.■  , /T"  U'  ♦ <?  /T"  )g- 

kq  qkq  qo  q q qo  kq 


*f  - «ufd-/Tdo,*d+‘cfkd/Tdo>*kdMcf/Tdo,*f+vf 


Herein,  v and  v are  the  d,q  components  of  the  armature  voltage  and 
d q 

v^.  represents  the  field  voltage.  The  time  constants  ‘1\  , l«d,q)  and 
the  dimensionless  c-constants  are  explicitly  expressed  in  standard 
parameters  and  listed  below. 


Ti  - V“V  - -V\v  °kd  ' -x,i/xd 


. lVx»’'1,d-xdl  _ lxd-x.i’  1>d-xdlx,i 
x.i  ix;,-x,r'x;;  (Xi-x,)2x- 


lxj-xd)xd  „ lxd-xt)(xd-xd)xd 
"dkd  - (XJ-Xt)X3  • °df  ' (Xa-X,)(X’-XtIX" 


<v*d,x* 


(X*-X"MX.-X  )X„ 
d d d £ 


ckdf  - (Xd-VX»  ' Cfkd=  2 

at  a 


(Xd"X£)(VX£) 


C * ( V 1 - Y \/Y"  . C = 

ckdd  1 d AC,/Ad  ' fd  (X'-X„)X" 

d t d 


c = -X  /X"  , c , = /(X  /X")'  - *X  /X" 

q q q q^q  q q q q 


It  is  noted,  that  an  additional  modified  armature  leakage  reac- 
tance X ^ has  been  added  to  the  conventional  set  of  standard  parameters, 
since  the  conventional  set  alone  does  not  completely  describe  the  machine 
behavior. However,  if  the  "leakage"  inductance  of  damper  winding  in 


the  d-axis  is  negligible  with  respect  to  its  self  inductance,  it  can  be 
readily  shown  that  X^  = X'^. 

A sudden  short-circuit  is  taken  from  a steady-state  no-load 
condition,  where  the  field  current  i^  is  adjusted  to  obtain  some  desired 
value  of  the  no-load  armature  voltage  E.  Thus,  before  the  short-circuit. 


0,  v = E,  v = v 

q f f 


and  if  the  elements  of  the  5-state  vector  i/<  are  ordered  in  the  same 


sequence  as  in  (73),  the  constant  state  is  given  by 


I 


with 


\J>*  <to>  =>  0 . (80) 

If  the  three-phase  short-circuit  current  i . , and  thus  its  d, 

abc 

q components  i^  and  i by  applying  Park's  transformation,  and  the  field 
current  change  i*  ( = i^.  - i°)  are  considered  as  the  observable  output, 
then  the  generally  noise  corrupted  output  equations  are  described  by 


id  = («/Xd)  (cd**  ♦ cdkd^d  + cdfi|)*)  + 4d 


1 = 


(u)/X  ) (c  ip*  + c . ) ♦ £ 

q q q qkqykq'  'q 


l'1/RfTio>  <cfd*'d  4 WiU  * VI’ 


+ c 


(81) 


where  K is  the  field  winding  resistance,  and  £.,  t,  and  £,  represent  the 
r d q f 

corresponding  measurement  noises. 

b.3.2  Machine  Parameter  Estimator 

The  standard  parameters  to  be  estimated  are  X.,  X',  X'| , X,.  X , X",  K 

d d u ( q q 

T*  , T”  and  T"  . Note  that  in  (81)  Is  also  considered  not 
do  do  qo  f 

known,  therefore,  it  has  to  be  estimated  together  with  the  standard 
parameters.  Equations  (7*0,  (80)  and  (81)  are  now  con- 
sidered to  correspond  to  the  general  formulations  in  (5S)  and 
(bb)  . In  (78),  -E  Is  regarded  as  an  input  whose  absolute 
per-unit  value  is  equal  to  the  no-load, per-unit  voltage  of  the  generator 
before  short-circuit,  thus  known  from  the  short-circuit  test.  The 
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remaining  equations  needed  to  construct  the  parameter  estimator  are 
obtained  by  applying  (58),  (59),  (62),  (64),  and  (65),  where  the 
expressions  for  3A(a)/3a,  3B(a)/3a  and  3C(a)/3a  have  to  be  evaluated. 

A block  diagram  for  the  estimation  algorithm  is  given  in  Figure 
13,  where  the  real  output  currents  i!jf,c(A)  and  if  (A)  are  expressed  in 
amperes.  Therefore,  these  currents  are  divided  by  their  base  values 
iB  and  ipg,  respectively,  before  they  are  compared  with  the  correspond- 
ing per-unit  currents  of  the  model  reference.  The  base  current  i„  is 

D 

easily  computed  from  the  machine  rating.  Depending  on  the  choice  of 

35 

per-unit  system  for  the  rotor  quantities,  a base  current  ifR  can  be 
obtained,  e.g.,  that  value  of  the  field  current  which  produces  rated  no- 
load  voltage.  It  is  noted,  however,  that  ipg  only  affects  the  per-unit 
value  of  Rp,  but  none  of  the  other  parameters  to  be  identified. 

6. 4 Simulation  Results 

A digital  computer  program  listed  in  Appendix  E has  been  developed 

for  the  WLS  and  ML  estimation  and  applied  to  a 120  kVA,  208  V,  400  Hz 

aircraft  generator  which  was  simulated  as  a fifth  order  model  on  a 

digital  computer.  Its  per-unit  parameter  values  with  a time  base  of 

1/2513  s were  taken  equal  to  the  nominal  values  listed  by  the  manufacturer 

as  Xd  = 2.10,  = 0.216,  X|j  =»  0.186,  X^  = 0.04,  Xq  = 0.786,  X"  « 0.105, 

R - 0.0189,  T’  * 522,  T"  «=  18.2,  T"  - 115. 
a do  do  qo 

The  short-circuit  takes  place  at  unity  no-load  voltage  (E»l), 
which  corresponds  with  a field  voltage  Vp  of  0.00210  and  current  i°  of 
0.0510  for  the  simulated  generator.  The  short-circuit  current  i^|u.  and 
field  current  ip  in  per-unit  are  displayed  in  Figure  14. 
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Figure  1A.  Oscillogram 

8< 


i 


The  parameter  estimation  algorithm  was  performed  digitally  with 

a sampling  time  of  0.2  pu,  whereas  the  total  number  of  samples  is  900 

corresponding  to  an  observation  time  length  of  180  pu.  The  gain  element 

gii  (k)  in  (1.2-12)  was  selected  equal  to  unity  if  ] Aou/a^  ) < 0.25, 

~k  ar  Ak-1 

whereas  g^.-  (k. ) equals  0.25cu  / 1 Aoc  ^ | if  lAa^/ou  | £ 0.25  for  all 

parameter  a^. 

1.4.1  WLS  Estimation 

Here,  the  weighting  matrix  W was  chosen  equal  to  the  identity 
matrix.  Several  sets  of  initial  parameter  estimates  were  subsequently 
tried  with  values  equal  to  40%,  60%,  80%,  120%,  160%  and  200%  of  the  true 
values,  and  finally  with  values  chosen  at  random  between  40%  and  200%  of 
the  true  values.  The  results  are  summarized  in  Table  4.  The 
behavior  of  successive  estimates  is  displayed  in  Figure  15,  which 
is  typical  for  the  entire  set  of  parameters,  except  , whose  behavior  is 
visualized  in  Figure  16.  Figure  17  shows  the  behavior  for  all 
parameters  with  their  initial  estimates  chosen  at  random. 

The  results  show  that  the  estimates  for  all  parameters,  except 
, converge  to  the  true  values  in  a few  iterations  with  less  than  1% 
error,  while  it  takes  some  more  iterations  for  X^  to  converge.  It  is 
also  observed,  that  the  convergence  is  faster  for  initial  estimates 
with  smaller  deviations  from  the  true  values. 

1.4.2  Output  Noise  Effect 

To  investigate  the  effect  of  output  noise  on  the  WLS  estimation, 
a zero-mean  gaussian  white  noise  with  constant  standard  deviation  of  5% 
of  the  steady-state  magnitude  of  ifl^c  and  i^.  was  properly  added  to  the 
output  of  the  simulated  generator.  Keeping  all  other  conditions  the 
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Figure  17.  Ratios  of  WLS  Estimate  to  True 
Value  for  Initial  Estimates  at 
Random 


same  as  before,  the  machine  parameters  were  estimated  using  randomly 
chosen  initial  values.  The  numerical  results  are  summarized  in  Table  5 


pw"",p’ 


and  displayed  in  Figure  18. 

From  the  results,  it  is  observed  that  the  convergence  rate  is 

about  the  same  as  in  the  deterministic  case.  The  estimates  for  parameters 

Xd>  X^,  Xj,  X|j,  Rg  and  T^q  converge  to  the  true  values  employed  in  the 

generator  simulation  in  a few  iterations  with  less  than  1%  error,  whereas 

the  estimates  for  X. , X , T"  and  T"  converge  to  biased  values  resulting 
x.  q do  qo 

in  errors  between  2%  and  13%  with  X^  showing  the  largest  error. 

6.4.3  ML  Estimation 

The  maximum  likelihood  estimator  was  then  applied  with  the  same 
conditions  and  output  noise,  whose  inverse  covariance  matrix  R ^ = diag 
(0.02205,  3.845,  0.02205  x 10^.  The  results  for  various  initial  estimates 
are  listed  in  Table  6.  For  easy  comparison,  the  convergence  behavior 
of  successive  estimates  for  randomly  chosen  initial  values  is  depicted 
in  Figure  19. 

The  results  show  that  compared  with  the  WLS  estimation,  the  ML 

estimator  has  about  the  same  convergence  rate.  However,  the  accuracy 

is  significantly  greater,  viz.,  all  parameter  estimates  converge  to  the 

true  values  in  a few  iterations  with  less  than  1%  error,  except  X^,  T^jo 

and  T"  , which  have  an  error  between  1.5%  and  2.5%. 
qo 

6.4.4  Input  Noise  Effect 

The  effect  of  input  noise  on  the  estimation  results  was  investi- 
gated by  adding  an  input  noise  to  the  generator  simulation  while  keeping 
all  other  conditions  including  output  noise  the  same  as  before.  The  input 
noise  was  generated  as  a zero-mean  gauss ian  white  process  with 


L 


85 


OinvOOlO'OOyDOyOO'OyOOON'OOin'O 

o oo4r^O\^>r^or^or-or-r>.ooor^f^or'r- 

0'00>0\v£)C7»C7>r^O>ooa>>OC7>C7>,^K)CJt(7»’90>(7» 

fNjOOr-HOOr-HOOOOOOOOOOOOO 


OromOOOinOmOinO**  inOrO*-HtOO>0in 
O OOoOOO^OOOooOaOOvOooOvnt^oOOr-HOO 
Z T)  OC7»(7»\£)a»a>fN*0>fl0CnyD0>0J^Nrv'0»i/)0K7l 

‘H  

CNlOOr-HOO^OOOOOOOOOOOOO 


Or-^lOOTjLOOlOOinOLOmOr- lO  O 0O  l/) 
0f^00C7»000000r-«00a>^00r-»0 
ocnovoa^orMOcx>osDoo^ta»oor^oo 

fN  O H fH  O H r-H  fH  O r-H  O rH  r—i  o O rH  rH  CD  *H  H 


CD  O *H  O O iH  CD  *H  CD  r-1  CD  rH  r-H  CD  fN  rH  rH  CD  rH  rH 

oaooooooooooooooooooo 

OC7>Ov0OOrMO00OvC>OOTTCy»OOOOO 


(N|  O .-H  r-H  rH  H H rH 


lo  r-H  r-H  j o O r-H  H (N  r-H  r-H 


c o-^cnor'Cr»ocr»ocr»oc7vcnocj>cr»cnocT»cr» 

o :croaoa>oaia>oo»oaioa>a»oa>wchoc7>a» 

H < X OC7>a»v0<^C7»NQ>00C7>v00»<7»,«tC^0»a>'0C7>0> 

ti  fMOO^HOOr-HOOOOOOOOOOOOO 


OrMOOrOOOOOOOOOOOOOOOiO 
cr  o ^ oo  o r-  oo  lo  ooooooooooo^qoooonoo 

0<  O <7>  C7>  C7>  a>)rs|  O)00CT)>£)O>C7i^NC^O^00a»(^ 

PM  O O r-H  o o|r-H  ooooooooooooo 


6 

•H  cH 

H-*  < X 

C/l 

UJ 


o ^ aiooaiooiooioc^o^oco^c^oioai 

O^OOOOOOOOOLOOOir>COOOrHO 

OOOOvDCTxOr-JOcXiOvOsOO-^rHTtOooOO 

fNO-Hr-HOr-HrHr-HOrHOOrHOOOrHr-Hr-Hr-H 


4):tj 

6 «* 


OTJ  OOaOOOOOOOOOOC7>rHOOOO 

oaooocr^ooooooooooooooo 

Oa>Os£>C7>OrMOC0OvDOOTtOOO'0OO 

CN  CD  r-H  r— H CD  r-H  rH  »“H  CD  r-H  CD  r-H  r-H  CD  r-H  r-H  r-H  rH  r-H  r-H 


OOoOOtT  cxDOoOOOOOrtOOOr-hOOOOrHOO 
or-  £7ioaio>o(7»oa>oa>0)OK)ooo)00»a> 
oo»ot\ocTka»(NO»coo>\oo>o>Tt(j)a>(^^o>c^ 

rjOOr-HOO^HOOOOOOOOOOr-HOO 


Or^NOO»f'JO(NOf^O'OMOCOO>NONM 
-a  OOOOCOOOOOOOOOOQOrHOOOO 
<x  OC7>Ov00>OrMOa0O'0OO'*fr0>OO(MOO 

rv|Or-Hr-HO*-HrHrHOr-HO»-HrHOOr-Hr-Hr-Hf-Hr-H 


I o m mo  m mo  tt  o o in  00(0  in  o o m oo  i 


88 


standard  deviation  of  1%  of  the  steady-state  value  of  the  corresponding 
state  variables.  An  identity  coef  f ic tent  matrix  for  the  input  noise  was 
used.  The  parameters  were  estimated  with  randomly  chosen  initial  values. 
The  numerical  results  are  given  in  Table  7 and  the  convergence 
characteristics  of  successive  estimates  are  plotted  in  Figures  20  and 
21. 

The  results  indicate  that  the  convergence  rate  is  practically  not 

affected.  The  WLS  estimator  yields  for  X.,  X„ , X , T'  , T"  , T"  and 

d ? q do  do  do 

T^q  an  estimation  error  ranging  from  2%  to  50*.  The  estimation  accuracy 
of  the  ML  estimator,  although  somewhat  degraded,  is  still  satisfactory, 
viz.,  the  errors  are  less  than  2.5*  except  for  X^, , whose  error  is  14*. 
Thus,  the  results  support  the  superiority  of  the  Ml.  estimator  over  the 
WLS  estimator  if  system  noises  are  of  significance. 

1.4.5  Effect  of  Short-Circuit  Resistance 

Since  in  real  world  implementation,  the  short-circuit  connections 
have  a finite  resistance,  it  is  desirable  to  investigate  the  effect  of 
a resistive  load  on  the  estimation  results.  Using  the  wu ighted- least- 
squares  estimator,  the  machine  parameters  wore  identified  with  random 
initial  estimates  for  several  per-unit  values  of  the  load  resistance. 

The  simulation  results  listed  in  Table  H show  that  the  estimation 
accuracy  is  affected  by  the  load  resistance  value  with  the  greatest 
accuracy  occuring  at  R^  » 0.0.  However,  a load  resistance  of  less  than 
Rj  - 0.5  pu  still  provides  the  same  accuracy  but  at  the  expense  of  a 
few  more  iterations.  This  is  very  favorable  since  in  practice  the 
short-circuit  connections  will  have  a much  smaller  value  than  0.5  pu. 


'»() 


Estiraat 


Number  of  ltei.it  ions 


Figure  21.  Ratio  of  MI.  Estimate  to  True  Value  for  Random 
Initial  Estimates:  Input  Noise  Effect 


6.5 


encountered  for  Initial  estimates  within  the  interval  from  402  to  200% 

of  the  true  values  employed  in  the  generator  model  simulation.  With 

the  same  weighting  and  gain  matrices,  similar  results  have  been  obtained 

for  different  sets  of  parameter  values  typical  for  large  60  Hz  machines. 

For  initial  guesses  outside  the  above  mentioned  interval,  the  algorithm 

converges  slower  and  may  even  diverge.  This  is  not  surprising  since  the 

estimator  equations  are  highly  nonlinear  in  the  standard  parameters.  In 

the  practical  implementation,  however,  initial  estimates  can  be  taken 

36 

equal  to  the  values  obtained  from  a reference  table,  design  data  or 
conventional  tests,  so  that  convergence  problems  should  not  occur. 

The  simulation  results  show  that  the  estimators  are  very  accurate, 
viz.,  in  the  absence  of  environmental  disturbances  and  with  the  availa- 
bility of  accurate  measuring  devices,  the  estimation  error  for  all  para- 
meters is  expected  to  be  less  than  0.1%  of  the  true  value.  For  noise 
levels  which  can  be  reasonably  expected  under  experimental  conditions, 
reasonably  accurate  parameter  estimation  can  be  anticipated  by  employing 
the  maximum  likelihood  estimator  which  utilizes  the  a-priori  knowledge 
of  measurement  noise. 

The  estimators  are  very  effective  considering  that  all  parameters 
are  simultaneously  obtained  from  the  data  of  a single  test,  where  it 
is  not  even  necessary  to  reach  the  steady-state  condition.  Note  also 
that  the  q-axis  parameters  are  extremely  difficult  to  determine  accurately 
by  conventional  test  procedures. ^ 
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\ 

The  study  on  the  effect  of  resistive  load  on  parameter  estimation 
indicates  that  the  best  results  are  obtained  with  a zero  load  resistance. 
It  is  also  noted,  that  the  sudden  short-circuit  test  provides  superior 
results  if  conpared  with  a step- function  field-voltage  test  with  short- 
circuited  armature  winding.  In  the  latter  case,  the  subtransient  behavior 
of  the  armature  currents  is  practically  not  noticeable  so  that  larger 
estimation  error,  slower  convergence  rate  or  even  divergence  can  be 
expected  as  indeed  has  been  found.  Consequently,  the  sudden  short- 
circuit  test  is  an  extremely  proper  choice  for  parameter  estimation. 
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APPENDIX  A 


l.OWKK  ROUND  ON  ERROR  COVARIANCE 


The  Inverse  of  Fisher's  Information  matrix  provides  a lower 
bound  on  the  error  covariance  of  the  maximum  likelihood  estimate 
(35,32).  A similar  approach  results  in  a lower  bound  on  the  error 
covariance  of  the  maximum  a posteriori  estimate  134).  These  bounds 
are  genera  l izat ions  of  the  Cramer-Rao  bound.  This  appendix  outlines 
derivation  of  this  bound  for  the  estimator  used  in  this  research. 

If  the  estimator  were  to  converge  exactly  to  a parameter 
estimate  y,  then 

ki  dh(t.)T  _ 

} l — ~r  V l « C t 1 - hfx.v.t  )]}  - V (v-m  ) - 0 . (A-l) 

, , dv  ~ w k k v v 

k-1  y 1 } 

In  words,  the  model  output  hfx.y.t^l  would  exactly  equal  the  observed 
output  This  condition  is  the  theoretical  best  estimate  of  the 

parameter  vector,  but  is  never  achieved  in  t he  presence  of  measurement 
noise.  It  is  reasonable  to  use  this  theoretical  performance  limit  tv' 
obtain  a lowet  bound  on  the  error  covariance.  Linearizing  hlx.v.t^) 
about  the  true  parameter  vector,  v,  gives 


V 

k>l 


dhit.  ) 

< I 

dv 


w 


h l x , v , t ^.1 


dh  i 

vlv  ' 


iv-vl  )) 


(A- 21 


-V_liy-yl 


V ' (v-ml 
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0 


')  ? 


_ — - 


Solving  for  (y-y)  yields 

k 


(y-y) 


, f dh(t.  ) 
-1 , r k 


R l 

k-i 


dy 


\>k>> 


(A- 3) 


where 


r-  a 


kf  dh(tk>y  _j  dh(tk) 


k-1 


dy 


w dy 


) + V'1) 

y 


and 


(A-4 ) 


wk  “ Z^lk^  ~ h(x,y,tk) 


Squaring  and  taking  expected  value  yields  the  desired  error  covariance. 
Equation  (7.1-3). 


cov(y-y)  - R \ l 
k-1  j-1 


k ...  ,T 

f f dh(t,  ) 


dy 


-w  d^h) 


V E(w.w . )V  , 

w k j w dy 


kf  dh(t  )T  ^ . _x 

-2  Y ~\  V ^[w,  (y-m  )Mv  + v Mr  , 

M dy  |~  w 1 kv  y y y 
k-i  i y 


Now  the  noise  samples  are  assumed  independent 


E(v, 


r v 

" 


V if  k-j 


if  ki*j 
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(A-5) 


(A-6) 


A 


T 

and  E[w^(y-jn)  ] is  assumed  to  be  zero,  i.e.  the  noise  is  assumed  uncor- 
related with  the  random  parameter  vector.  As  a result,  the  desired  bound 
is  simply 


cov(y-y)  * R *RR  * R 


(A-7) 
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APPENDIX  B 


NOMINAL  PARAMETERS  FOR  SIMULATION 


The  nominal  parameters  used  in  the  simulated  experiment  of 
Section  II were  approximations  derived  from  nameplate  data,  steady- 
state  measurements,  and  a list  of  typical  machine  parameters  [3]. 

From  the  generator  nameplate  the  rated  values  of  armature 
voltage  and  current  were 


= 132.79  volts,  line  to  neutral 


(B-l) 


*B  ‘ I3§7§9  ‘ 7'531  amps  ' 


(B-2) 


The  angular  speed  of  the  equivalent  two-pole  machine  is 


id  - 1200  x 77  x 2 = 251.33  radians/second 
o o u 


(B-3) 


The  results  of  steady-state  open-circuit  and  short-circuit  tests 
are  plotted  in  Figure  B-l.  The  rated  value  of  open-circuit  armature 
voltage  corresponds  to 
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. 31  amps 


(B-4) 


1 


FB 


neglecting  saturation.  Thus, 


FB 


VB1B 

lFB 


3226  volts 


(B-5) 


and 


~ 4.87  x 10  per  unit 


(B-6) 


Rated  short-circuit  current  corresponds  to 


1 - .253  amps 

r d 


( B—  7 ) 


therefore , 


XFS  1_ 

I„D  w 

IB  o 


3.25  x 10  per  unit 


(B-8) 


A typical  value  of  transient  reactance  is 


X’  - .322  X.  . 
d d 


(B-9) 


Therefore,  from  the  definition  of  transient  reactance. 


X - x-  - 
Xd  Xd  Xf  * 


(B-10) 


and 
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10. 78  x 


(B-l  1 ) 


‘df 


Lf  - 678L 


10  3 per  unit 


The  remaining  d-axls  parameters  are  roughly  approximated  by  assuming 
equal  coupling  from  stator  to  all  rotor  circuits. 


Lkd  “ k"Ld  “ .<>78x3.25x10  3 - 2.20xl0'3  p.u. 

Lfkd  “ ^^fHd  " ^ ^df  ” kLdf  “ 4 . 01x10  3 p.u. 
Typical  q-axis  parameters  give 

L ~ .6521.,  - 2.12x10  3 per  unit 
q ™ d r 


(H-12) 


(B-l 3) 


CB-14) 


and 


■ L - L”  - . 55x2. 12xlO~3  « 1.17x10" 3 per  unit  (B-15) 


With  the  rotor  at  standstill,  direct  current  measurements  give 


Rn  - .2511  and  Rp  - 24012  , (B-16) 


or 


Rj  ■ 1.418x10  ^ p.u.  and  R^.  - 2.30x10  ^ p.u. 


(B-17) 


The  damper  resistances  are  not  measurable;  however,  typical  time  constants 
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467  sec. 


T ' - — 
do  Rj 


and 


T.  " ~ .02  T.  ' 
do  — do 


The  result  is 


*kd 


t fkd 

hcd  " Lf 

..  <r  M 

do 


.00935  sec. 


= .0759  p.u 


For  lack  of  better  information,  we  assume 
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APPENDIX  C 

PARAMETER  IDENTIFICATION  PROGRAM 
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SU  9-0  JUNE' 3 USED  XSS7,  G-A07, 
CVe.ii.AT  LVlI  HAS  XES7.IKV  AND 
CoMPILlR  F-'Et  . DOUPL.-  PA  f CliXON 
COMMON  EF  < , RL  K » WR 


PROG - A M 


F.RCM  o Aw  A MET  EKC 

S ■ J l V 5 * DCHY7,  DF7  , AF7,  UDAT  7,  INV 
30uV5!  OVERLAY  0VL2  OCH Y7 , UO AT  7 , GR AO7. DF7 


DIMENSION  Y(11),ZF(4),'(5),X1(F)»W1(11,4),Q(4»c),W2(11»11) 
DIMENSION  M(ll  ,111  ,WM 11) ,C0 (111 ,C  (4,11) ,X2<5) .MM (11, 12) 
DIMENSION  OX  < ? » ,0X1(5)  ,DX2  (5)  ,W3U>  ,Z(4)  ,W4(11> ,0M<4, 11) 
DIMENSION  OXY (5, 11) ,0(4,5) ,OH (4,11) ,OE  Lllll.DY  (5.11) 
DIMENSION  PXY(F,11),PX0(5,5I  ,P(11)  ,VM(11) 

C cXTE-NAL  DV-1.0VL2 

C THE  FOLLOWING  IS  EOUIVAcENT  TO  THIS  EXTERNAL  STATEMENT 
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OVwl=0 
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-.EADU)  EFK.Kl  1,K1,SL2.K2,  K3 
JO  Z |_=1.NZ 

7 Rt uO(-)  (Q<w,M).M=l ,NZ> 

•vEmD(w)  (P(L).L-l.N) 

- EA  0 (*- ) (Y(J),J=1,N) 

K = 1 

DO  A K=1.KF 

F;AO(E ,1002)  (Z(L) ,L=1 ,NZ) 

WRIT*  (7,200  51  K,  (Z <Lt ,L  = 1,NZ) 

A WRITE (fc, 2003)  K,  (Z(L).L  = 1.NZ) 

WAIT1 ( fc  » 2 00  a ) (Y(L».L=1.N) 

JO  IN  J»1,N 

19  Y.N(J)  = Y(J)  __ 

SOV-O  ” 1 rnfmmrtmmt 

1 = 0 
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WN(L) =0. 

LDAO  F I -\ j T OVERLAY  CVLl 
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C AwL  X Si7 (X.NX.Y.K.CO.LY.Th) 
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COM  Me  NT  1 


COMMENT  I 
CO 

COMMt NT  I 
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SO 


90 


PXY  ( J * L > = C Y ( J » L ) 

00  12  L=1,NX 
00  12  J=1,KX 
PXO ( J»  L) =0  . 

00  13  J=1.NX 
PXO  < J,  J) =1. 

RlWINO  5 
T<=0  . 

<=1 

LOAD  7 i Nr  LOOP  0 Vt  PLAY  * VL2 
C A„L  OV. JO< 9.CVL2, IOV, IGMl 
IF  (IOF.U.1)  STOD  "iPFO*'  IN  LOAD  DVl2" 

OEGIN  TIME  LOC’p 
CunTInUI 

StJ  Jp  LIU-TICMSt  0 = [H-/DX,  C=OH/0Y,  OXY=OX/JY 

IFCK.EQ.Kl)  F L <=c  L 2 

IFCK.cQ.KZ)  rL<  = -Ll 

IF(<.i_0.<3)  L < = p L 2 

TK=T« (<-l) 

JK  = 1 
TJ=T/NZ 

IF< J<. GT.N7)  GO  TO  93 

CAwL  DC9Y7<X,NY.Y,fc»l7.CO.C.O.Z.T<tTH) 

00  AO  i.=  l,NX 

Oo  SO  M=l,v 

OXY  (L,  M)=PX  Y (L.M) 

00  30  J=1»N,< 

OXY  CL  » M)  =CXY  (L  .Ml  ♦rxo  ( L . J)  * 3Y  ( J,  M ) 

00  90  L=1,!.Z 
DO  90  M=l.*>. 

DMCL.M)=0(t.M) 

00  30  J=1,‘X 

OM(L.M)suM(L.H)»D(L.J)*0>YCJ»Ml 
0 AwL  UC  AT7  ( X ♦ K X . Y ♦ l.«  PX  Y ♦ FV  0 . 7 J) 

ZP<  JO  =Z  ( JO 
DO  31  Ms  1 , N 
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■■IMP 


91 

93 

100 

120 

150 

170 

190 


OH( JK,h )=OM(JK,M> 
TK=TK+T J 
JK= JK* 1 
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GO  TJ  92 
COnT INUE 

WRITE <7,20031  K, <ZP(J» ,J=1.NZ> 
DO  1U0  L = 1,H 
DO  100  M=1,N7 
WI(l.M>=0. 

JO  100  J=1,NZ 

W1<l,M)=W1(l,R>*OH(J,.L>*J<J,R> 
uO  120  1=1,9 
DU  120  M = 1,N 
W 2 ( L , M » = 0 , 

JO  120  J =1  , NZ 

W2<L,M)=W2<l,M)+W1(L.J)»UH<J,M> 
00  150  „=1.N 
CO  15  0 M=1,N 
WM(L,M  >=W3(l,K  >*W2  <1 ,M) 

READ  (5 ,1002)  < Z ( L » tl * 1 ,NZI 

DO  170  L=1,9Z 
W3<L)=Z(u)-Ze<U 
OC  190  „=1.N 
W-»(L)  =0, 

JO  1 9fl  J = 1 * N 7 
W-»(_)=W4(L)  ♦ W 1 (L*J>*W3(JI 
JO  200  »=!,.< 


200  WN(w)=WN<l>*WMU 

WR*Tf  ( 10)  K,'’<13>M 

<=<M 

:f(k.lc.<f)  gc  to  ao 

COMNl  NT i END  OF  TIM:  LOOP 
COMMENT  I LoAD  LVLl  AGAIN 

CAlL  OVLJJ(3»CVL1.IOV.:Ok> 

IF  ( 1 OF  • 9l • 1 ) DT^  "fRf-0.\  IN  LOAD  OVtl” 

c 

OO  699  _ = 1 . N 

699  MM(c,l)=WM(L.L)*c  (D 

UO  799  u*l,N 

799  WN(L»=WN(l)-R (L» *(Y<L)-rM(L ) > 

JO  600  J=1,N 
UO  610  „=1,N 

610  WM< J.L ) =WM < J ,L >*SC A 

W N ( J ) =WN( J) • SC  A 

600  WM(J,J)=WM(J,J)**SCP 
COMMcNTl  SJLVt  S i M'JuT  A NEol.'S  FQLATIONS 

DO  601  J=1.N 
DU  601  1*1, N 

601  W(J,_> =WM(J,l) 

OC  60?  J=1.N 

602  WM(J,N*1) *WN ( J ) 

Ms).*  1 

CA>.L  SOL Y5 ( WM » Cfc L » N, f ) 

CO  ?19  4=1, N 
W4<L) =0. 

UO  719  J=1,N 
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tin 


720 

I 


251 

COMMENT  t 


COKMf  NT l 
260 

262 

COMMENT « 


W 4 ( L ) = W 4 ( L ) ♦ w (l«J)*DEL  (J) 

WFITt  (F>  “WM(L.w)=“ 

WPuTf (6.20041  (W(L,L» .L=1.N> 

WRITu(fc)  •,0t.4  = " 

WPiTE (6.200*.)  (Of  L ( J ) » Js 1 » N) 

WR«.  n ( 6 ) "WM*  u i L = ** 

WiiTv  (6,2004)  (W4( J) , J=1,N> 

Wis.Tt (o)  "WN  = 

Tt  (6*2  00**}  <WN ( J)  « J- 1 *N) 

00  720  L=1,N 
Y (L>  =Y  (>)  ♦ DEL (L) 

WRITt  ( 5 ) "PAR AMLTtFS=" 

WrilTl  (6.2004)  (Y(l).L=1»N> 

MRiTf(lO)  “3AP AMLT[ ► 3=“ 

WRiT‘~(  10,2004)  ( Y (U  ,L  = 1 *M> 

CAlL  O V:  -»F  CJ  F > 

IF  (IGF.l'J.I)  c'AUSk  "Floating  point  OVtK FLOW* 
IF  (IJF,.i,il  PAUSE  "FLOATING  POINT  UNDERFLOW" 
IF  (I. LI. IT)  GO  t0  in 
LHJ  of  i te  :at i on  loop 

ACCEPT  ‘No.  CF  ADDITIONAL  I T i P.  AT  * JNS  = " , I A 
IFdA.LT.  1)  G r-  TJ  260 
: r = ir4  ia 
GO  TO  251 

jUlV l FjR  F I R A l VALUt  CF  X AND  ZP 
f LK  = R’L  1 

CAlL  XSS7(X,NX,Y,n,C0.DY.TH) 

7H  = TH*WR*  TO 
call  INV(W.W2.N) 

DO  2 62  »=1,N 
hZ(L»l) = W2  <L  «L  ) * 3C  A 

W R.  a T F ( b ) "ERROR  VAf-lLNC..  l ST1MATES" 

WRITE  (6.2004)  (W2(l.L> ,L=1.N) 

LlAJ  Oi/l 2 AoAIN 
CAlL  0VlJD<8.C VL2. iCV.IOR) 
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450 


4b0 


470 


4 BO 


410 

252 

250 
100? 
20  0 2 
2 0 0.1 
2 0 0 ■* 


*F(*  . 1 > _ Y 0 p " * F *■  i $ I ‘I  LOCO  UYL?" 

T J=T /HZ 

rfrvllr  <0,2  002  ) 

JO  2-0  X = 1 . <F 
IF  (K,r  Q,<U  -'L<si-L2 
IF  ( A . t_  Q.  <2 ) L<=^Ll 
IF(<.-Q,<J)  -\L  < = c L 2 

T<=T • (K-l ) 

DU  2 i?  J<=1,N7 

CALL  OCHY7(A.M*-,  Y.N.K7.C0.C.Dt7.  TK,TH) 

T <-T  K ♦ 1 J 

C A _ L G F A J 7 (AtNXyYfKfCX) 

do  -*5  0 j = i,n> 

A 1 ( J>  = CK  ( J)  • ; J/J.4-X  ( J) 

CALL  jcAJ',(XltNX«Y«NltDXl) 
oo  -*6  0 j = i «’4.< 

XI ( J» = Ox 1 ( J) ♦ DX ( J) » j/6. *x  ( J) 

CA*.l  3rvA37(Xl.KX,Y,N,0Al> 

DO  <*70  J = 1.N> 

X I (J)  = (OX  1 ( J>  * .1.  o>  ( J)  >*-.  J/8  ,+X  ( J) 
o A L « o'\l  J / l X l » HX  * Y » N , G * ? ) 

DO  4 10  J=l,  IX 

xl  ( J»  = (0x2  U)  • 4.  -DX1  ( j)*  S,  +jx  ( J)  ) *TJ/2.*X  ( J) 
CAwL  o A07  (M,  NX  ,Y,N,C>  1) 

JO  4*40  J = 1».A 

X(J)=(1)X2(J)*4.  + DX1(J)4  0X(J))*TJ/6«4A(J» 

ZP( JK) =Z < JX> 

H*IT. (7,2003)  <,  <7P(U  tL=l»NZ) 

W rX  1 T ’ (6,200.1)  K,  <ZP(l>  ,L=1  ,NZ) 
f o-  'i a t ( r.G  16 .e ) 

FOX^AT  (J<,,,\"lOXi"Z(M"l 
F Jr.  MAT  (lX»l4*l'Si  15,5) 

FO  <SAT  ( IX  ♦ IPS..  14.  ) 
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CON»:lE^  31 U °l  *'  PRECISION 
OViRUY  JVLl 

SUBROUTINE  XS£7(Y,KX» Y,N,CO,DYtTH) 

COMMENT  l CCMPUT  •_  ST- A OY  STATE  USING  CUR'-TNT  PAR  AML  HR  f ST IMATE  Y 
COMMENM  JY=-a  T.  JcK.  0 P >0  WFT  Y 
COMMON  cFK.RLKtM 

JIMS  NS  ION  X (*;>  > , Y (M  tOY(NX  ,N>  »C0  (N> 

11=1 
12  = 2 
x3  = 3 
I 4 = 4 

15  = 5 

16  = 6 
: 7=7 

16  = S 
:*= 9 

R = Y(U)+RlK 
S=oQ>.T  (2./J.) 

►>  = 3.  l*>15-}26535  697  1 

=>  = 2. **/3. 

OL  = Y (I 1»  *Y (16) 

QL=Y(I1»-Y(I6) 

OX  = Y <:«)  * (R*v2+w  “2*(?L*  Ql) 

X <:i)=tF<*YC2>*R**2/D> 

X (12  ) = ["FK*  Y ( 1 4 | / Y (19)  -YCZ»*«2*(JL*  ..FK*  W**2/0X 
X (13  )=r  F**Y  /Y  (T9) -Y(_2)*Y  (I3>*  OL'  EFK*  W«*2/DX 
X (.‘♦)=-cr<*Y  (I2)‘*QL*f:*W/D> 

X (15  I = -L  “<*  Y (-.2)  * Y (17)  * *•  * W / 3X 
DO  101  :=1,N 
00  101  J = 1 , N> 
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101 


COM'Il  --it 

160 


oy< j,: >=o. 

3YC  1,  ID  = -X  C 1)  *Y  (29)  • w*>*2*  ( U*Ll  >/0X 
OY < 1 1. I 2 > = » (1 1 >/Y ( 12) 

QYCl,  lb)  =-X  <1  1)  *Y  (T°>*  W**  2*  (IL-OU/JX 
cnii,;H*2*<ci)*(iA*y*Y(;<)  /ox  > 

DY(ii.jq»)  =-a  r i)  /y  (19) 

3YC2,;i)=k**4«Y<:2)**2*Y(H)*rFK*<Q*TjL)/CX**2 
- Y (I ?) • * 2*  LFK* W* • 2/OX 
2Y(I2.I2)s-2*V(;2)«dc*EE,<*W**2/0X 
OY{ = .FK/Y (I  9) 

JY  ( 12. 16)  «W**«.*Y  C 21  *♦  2*  Y C9l*c.FK*  (Qw-uul  /OX** 2 

♦ YC2MV*tFK*W»*i?/:x 

□ Y (12*26)  = 2*  Y (12)  **2*Qlv)-*  Y (19)*  LF<*H* *2/OX**2 
OYC2f.-J)s-XC2)/YC:Q) 

)Y ( 13,3 1)  = W**<-  •Y(:?)«V<:3)»Y(I3)».  FK-dL*  ( QL-DL  ) /OX* *2 
-Y(I?)*Y(I3)*.  FK«  W**2/CX 
3Y(I3,:2)  = -YC  3)  *2L*l  rK*W«  "2 /CX 
UYC  1,:J)=-YC2)*0l*EFK«w*  ‘2 /OX 

jrcj,i5)  =:‘f</y  C'i) 

0Y(:i,ib)=»Y(:?)*Y(I3)*j.PK*K-**2/3X 

♦ w ‘i’YUi’YCJint^i*'  r<*  jl*-  ( ic-ro  /ox**? 

OY  ( 1 1 , :«  ) =2*  Y ( :2)  * Y ( A 3 '*'0l*"*Y  (i-JI  *LFK*tf**?/DX**2 
3 Y ( I 3 . I 1 ) = - X ( 1 3 ) /Y  ( I 9 ) 

OY(iH.il)=-x  C -)  * W**2*<rjc*DL)*YC))/OX*X(I-)/UL 
o y ( : u , : 2 ) = . co/YC?) 

?Y(:-ti3)=-x(;‘*)*w**2*Y(:9)*(dL-0L)/JX-A<i‘.)/QL 
JVC  ♦,!■*»  = -X  <_  3)  ' < - 1/-  + 2*i  ■ Y (I  3) /OX) 

CYC-,.  >)  =>X  <:  4)  /Y  CO) 

[Y(c1:i)="(i:ri,^,?MQi.tjj*Y(:Mi/ox 
OY  ( * 5 , 12)  =>  Cc  ) / Y C?> 

OYCf  ,:6)=-a  c r ) *iw* *2*  ( c-'l)  *yc  n/^x 

i)V(ISi7|sH.r  ) / Y,(  : r> 

jy(if,:j.i  =-\<:c>  m-i /'  ♦?*■  -y  c>i)/o a) 

)Y(13.H)  = -X  C C I / Y (13) 

T H - :-C-IAu  XOTJk  ASbL-  C 0 = 0,>  l VhT  I V:  CF  TH  WkT  Y 

TH-A  ' AM4*U_/-  ) 

OJ  It')  M=1.N 
20(*<>  = 0. 

:0(Il)=.^*"/(<*'2*»**?*0w*,2) 

20  <:- ) =-20(11  ) 

20  ( I *)  =-r(*dL/  <-  *»2*W**  2*QL**  2) 

: jan 

;.NJ 
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r 


A<K=  A ( K.  t \ ) 

= 1 

00  to  Ki=KLl.*. 
ak*<=m  <k;  ik ) 

DC  10  L=<»l.N'l 


10 

A(<l,.)=ACiI,w>-<AK:K/tnK)*a(<tU 

D J 11  *1=1.1 

11 

Oc 1 A=A  C 1,; 1) 

C 0M'l?'<T  1 

s*a-':nd  ct  ftrrror,  scivr  fc  t o 

C(M  = A ( N » '.  ♦ 1 ) / A ( N » M 

Do  0 0 1=1. UNI 

0(  .-I*  = M >4-1  •»■•♦!» 

CN*=C(K-i) 

DO  JO  J = 1 « I* 

JO 

cn.  = c <:-o ( j*n- i»  * a (n-i  . j+f ) 

0 MI/A  (K-I  , >1-1  ) 

jci  «PrrA*-<i*i,iM) 

20 

uc-t  **p»r***»c’  ^ a 
mm-  u.ion  p-'ta 

101 

f ( i < , "jr: T t nant  = '*  , 1*1  ir.5> 

■»  £ 1 U * i i 

•_  .V) 
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i 


C&MPlLt’.  F-fcf  .DOUBLE  FflCISIUN 
OVr-luAY  OVl2 

j JbROJT  I NE  DCHY7  (X  , NX  , Y , N,  NZ  , C 0 , C,  D,  ZP  ,T,  TH  > 

COMMON  tFKt^LK.W 

OiME  NS  ION  X(NX  ) , Y(M  ,CO(N)  ,C(VZ,N)  » Z ->  ( 9Z ) , 0<NZ,NX)  ,A<5,5).F<5.5>,OT  H<4> 

‘jeal  io.iq 
:i=i 
12  = 2 

13  = 3 

14  = 4 

15  = 5 
I b = 6 
17  = 7 
I8=B 

3 = 3.1l»15926535'!J79 
a=2.*»/3. 

S = £TFr  ( 2 . / 3 . > 

OL  = Y (ID  *Y  C I fc » 

JE  = Y (I  1) -Y (16) 

CALc  AF7(Y.‘(.t.X.A.F) 

10=0. 

OC  40  J= 1 . 3 

40  :o=io+A(ii, j)»x(j) 

13=0. 

JO  50  J = 4,5 

50  I J=la*  AC4.  J)  » X ( J) 

00  90  1 = 1.1. 

Oj  90  J=1,\Z 
90  C(J.D=0. 

oo  9f  ; = i.t.x 
JO  95  J=l.f.Z 
95  0 ( J. I ) =0  . 

DO  100  1=1.3 

0(ll.II=3»A(Il,r)*CUS(w*74:H) 

0(*2.I)=3*>.(Il,:)*C0S(W*T  + TH-3) 

0(.3.;)=s»A(Xl,I)*CCS(W»T4rH4C«) 

100  u(-4,: ) =4 (12, I ) 

Ou  105  1=4,5 

o(ii,:)=-4"A(:4.i)*sif‘(w',T*TH) 

j (.  2. ; ) =-3-4  c*.,  1 1 *sif:  <w* ; +*h« 


J 
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in? 


no 


i?o 

130 

30 


uC3,.)=-.a*4(.4,l)*SiN<U"‘!*TH4t:) 

DO  110  1=1,;  U 

)=0. 

DO  110  J = 1,‘1X 
7PC>=Z0C)  ♦OC.  J>  *X  (J> 

00= Dl*  (V  (It ) *Y  (III  -Y  (I5)-**2)-Y<12>  * < V « 13»  *Y  (I2)-Y(i5)*Y  (I  3)1 
-Y(I3)*(Y(:h)»YC3)-Y(I2»*Y(I‘>I) 

OQ=QL  * Y ( 1 7 ) -Y ( 1 7 ) * *2 

01  = A(il,:i)*M:i)*ACl,I2l*A(I2)*A(Il,I3>*M:3) 

D4  = ACu,I.»)*X{rM+A(:4,iS>*K(If) 

C(il,i  1>  = -AC1  ,I1)*01 

C(11»A2)=(-Y(i  3)  *,<  (12)  ♦ Y (I  5 ) *X  (13)  l/DU-2*  A ( 1 1 , 1 2 » * 01 
0<;i,..3>  = <Y<:i>*x<ii>-Y<i2i*>(i2>  + Yii5>*x(:2»-Y<:4>*x<:3>»/DD 
- (A (13,13) +2- A (11,1 3 > ) *01 

CC1,I4)  = (Y(I3)*>(I1>-Y(:3)*X(I3)  )/0J-A(l2.I2)*ai 

C IJ  1 b)  = (-2  *Y  (I6)*X  (1 1 ) *Y  (I  3)*XC2)  *Y  (I2)*X(I  3)  ) /DO-2*  A (12, 1 3)  *D1 

c ci  ,;&i  =c  ( ii,  in 

UC4,I  1)=-A(l4,iu)*r<. 
c e»=-c(i  4,ii) 

C Cl  4, 1 7)  = C C4  )-X  (I?  ) ) /DQ-  ( A (I  5, 15)  *2*  A(14,I5>)  )*D4 
DC  120  1=1, V 

0 ( . 3 , 1 )=S*  (C(11,x)*DLS(W*t^TH*C)-C(14,  i)  • SI.'I(H*t*th+'3I  ) 
C(i2»I)s5MC(Il,II*CCS(W*7*7H«cl«C(I>,»  I)*SiN(W*T4  TH-O)  ) 

C Cl  ,1  > =S*  (CC  1 , i >*Ct  3(H*l+TH>-C(i*,D*SIN{W*T*TM)) 


DO  1 (0  i = l, N 

c(i«,,:)=o. 

D2=A(I?,:i)-X(Il)+A(I2,12)'X(12)+A(I2,I3)*X(I3) 

C(*4,:i)  = (YC3)*XU2)-Y(*6)-»>(T3))  /OJ- Ad  1,  ID  *02 

C(14, 12)  = (-  YC  3 )*X  (ID + Y(I  3)  *XC  3)  I/D0-2*A(I1,I2)*02 

C (i*,,*3)  = ( ( Y (1  ?)  -Y  (12)  )«  a(I1  ) * (OL-2*  Y ( i 3)  )*X(I2)  *Y  (I2)»XC3I  )/OD 

-(AC  3,13)  *2* A (11,1311*22 

CC-,14)  = -AC?.I2>«02 

CC4,I5)  = (Y  C 3 ) * > (Il)—2i_*>  (13)  ) / 00-2*  A (12,13)  *02 
DC  30  l=1,N7 
OTH ( l ) s 0 » 

0TH(  *1 ) = -S*  ( ID* SIN  ( W*t  *Th|  +IQ*COS  (H*  DU)  ) 

OTH ( 12) =-S* (10* ilM ( W*T+TM- PI +IQ*CO 5 (W*  T*3H-3) ) 

OTH ( : n = -s'  ( iu*si n(m-  T4-;m*»)  4io*oos(w*  r+i  hd)  > 

DO  35  1=1,3 

c(l,iii=o(l ,:i)*oth(l>*co(11) 

C (u  * 16  ) =C  C , If  )*DTH(L)*C0C6) 

C(w«It')=CC,I£>*D7H(L)*C0Cf) 

A f 1 U - I! 

wNJ 
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OOMaiLif  DCUBlr  o-lCiTJON 

b JF-JOUI  ML  0F7(DFX,DFY,X,Y,HX,N) 

COMMON  lFX.MLK.W 

DIM£nJ13-J  XI  MX)  ,y  (N)  ,CFX  (NX,  NX)  ,i>rri  (X  ,N>  « A (E ,r  ) 

Z 1 = 1 

T : y 

* U t 

i 3 = I 

1 4 = t 

15  = 5 

;fc=& 

: 7=7 
X a51) 

I9=* 

110=10 
A 1 1 = t 1 

o awl  AF7ir,.-i,Kx,Atrrx) 

DO  200  X = 1,M 
DO  300  J = 1 » N/ 

JfY ( J, X) =0  « 

> =r  c !)+-lk 

Ow  = Y (. II  ♦Y(Ic) 
lL*YCl)-Y(Io) 

0 = CL*<Y(:4|»Y(I3)-Y(:c)*»?>-YC2)*<Y(I3>*Y(I2)-Y(I5>*Y(IJ)) 
-Y(I3)MY(.4)*Y(MI-Y(12)*Y(I5)) 

jl*UF.\  {*.  1 , : 1)  • <(  , 1 ) +PF>  IT  1 . T2)«X  (12)  + UFX  I il  .1  3)*X  I)  3) 
0FY(u,u)=-Ani,ii)*:'i 

OFT  ( : 1 J)  =-  * ( Y <:  ?)*;  (12)*  Y C 1 *5  I *>  (13)  ) /0-2*  M 1 1.  13)  • 01 

JFY<;i.Al)=-.*(-Y(:M*X(;?)+Yir<*)*X(:l)-V(I*f)*XlIl)+Y{I7l»X(I3))/C 

-(A(M,I3)*2*A(Il,M))*Oi 

DFY(il,i4)=k*(-Y(I3)»X<:i)*Y(M)»X(M>>/0-M12.I2l*01 

OFY  (II  ,I5I=<*(2*Y(I5I*X(I1>-Y(IJ)*X(I2»-Y(;2>*YII3>> /D-2*  A 1 1 2 , I J)  • D1 

DFy(I1,I.))=D1/6 

- F = Y (10) 

J?=3FX (12, ill*X(ll> Ft FX (12, I2)*X(I 21 OFX ( 12, i J ) *X (13) 

OFr  112,11)  = KF*  l-Y  C 3>*M.2>  ♦YC5)  *x(i  J ) ) / 0- A 1 1 1 , 1 1 • *02 

l)FY  (12,12)  =KF*  (YC3)*>(Il)-Y(M)*X(U)»/0-2*A(Il,I2»»02 

OFY  ( 12,  13)  =-.F*  ((-Y  (I‘.  ) *Y  C 21  )*X(I  II  ♦(2*YCJI-DL)  *X  (I2)-V(I2)*X  113)  >/D 

- (~  I J.(,:3)  +2*A  (il  ,!3  ) ) • 07 

DFY(I2,I‘*)=”A(*2, 121*02 

OFY (12 ,151  = \F*  (-Ydlincil  ♦OL'X  ( MII/0*7*i(I7,I  I >*03 
OFY (13. I ») =02/  F 

• XO  = YC10) 

D3=3FX  (I  1 ) • X(  .1  ) +[jF>  (;.(,  131  *X  C3I  ♦LUX  ( 1.1,1  JI*X  C 1) 

JFY  ( M,:  1 > =VKfJ*  (Y  ( Jc>  •>  (12)-YIT-I  • Ml  J II  /0-A(I  1,1  11*0  J 

OFYC3,l7)='-Xu»(-Y(Ic.)*X(;ll-Y(:(l*X(l2)»7*V(l2l*x(:JI)/T-7*A(Il,I?l*03 

ofy  ( i i ,i  n - -c  k n - (Y(:u)*x(:i!-Y(;2i*x(i2M/o-(fl(r3.M)4’*«(it,i3M*OJ 

jFY  m,I4»  S^KC  * (Y  (i3l*X(Jl  )-rL*X  (Ml  l/0-«  (I7.I?)*C3 
OFY  ( 13 ,15)  = <KT  * (-Y  II?)  * X (i  l ) ♦0L*X  (i?  I ) /O-  3*  A ( 12,  I 3 MU  3 
OFY  (11,110)  = 03A:K0 
0 1=(}U  * Y ( . 7 ) “Y  I : 7 ) * * 7 

0“rDFx  (*■*,:<*)•  XCI. ) »DF  Mi -..IE)*  X ( l'>» 

ofy  ( i- . : *5 » =a  t j :<*)*r4 

OFY  (I*,  ,i/l=>»*(X(I5)«XlT4)  ) /Jll-IA  (1E,.5  A(I4,IF  ) ) * 04 

0FY(l4,Ii»=D4/“’ 

<XQ=Y ( lit) 

os =ofx  (ij,:9i*x(:4)*i.Fx(.‘,i5»*x(;5> 
ofy  ( 15,ifc)  = <KO“X  (151 /00*A(  ,4,14)  •lV-, 
ofy  (15 ,17)  = kkq*x  (i  «•)  /o'i-  (A  i ,:m  4i  14,  :r  i i*o5 
ofy  ns,;  in  =co/‘»k.i 
w ; u ? m 
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t i^jU  COPY  FURNISHED  TO  DDC 

COMPILE  OJ'JOLE  PPt  CIS  ION 
SUBROUTINE.  AF7  (Y.N.NX.A.F) 

COMMENT » THIS  SUSROUTm  COMPUTES  A=IN(ML>  AND  F=-R*A 
COMMON  EFK,4LK,W 

D I Mi  NS  13 N Y(N)  ,A (NX  t NX  > « F ( NX  » NX)  ,*(=> 

11  = 1 
12  = 2 
1 3 = 3 

I 4=  4 

15  = 5 
ib  = 6 
17=7 
18  = 8 
10=9 
I 1 0 = 1 11 
111=11 

DL  = Y (1  It  +Y (lb) 

QL  = Y Cll-rilbl 
O'j  ?o  : = i.nx 
JO  30  J= 1 « NX 
30  A(J,:t=0. 

do  4tt  : = i, nx 

Jo  4 0 J=1,M 
40  F(J,1)=0. 

D1  = Y (1 4) *Y (I3> -Y  (I5)*»  2 
J2  = Y (i 3)  *Y (I  2) -Y (15) *Y (13) 
l/3  = Y (i4)*YI*3)-Y(;2>nCS) 

0 = DL*U1-Y  C2)*02-Y  (13)  “03 
DQ=Ul‘*Y(17)-Y(I7)**2 
A ( j.  1 « 1 1)  = 01/0 

aci,:2»  =»l.2/c: 

A (1 1 ,13)  =-t  W' 

A (12  , : 1)  =A  (11,12) 

A (*2t  - 2)  = ( DL  * Y (I  3)  - Y (I  3)#*  2)  /O 
AC2,:3)=-(JL*Y(*5)-Y(I?l»V(13))/j 
AC  3,1  1)  = A ( 1 1 » ? 3 ) 

A C 3 * * 2 1 =t (12. 1 3) 
A(l3,:3)=Cu<*V(I«)-Y(I?)^*2l/j 
M.4iM=yC7l  /Df3 
A(_~,lr>  s-YCD/DQ 
A(I5*.h)  = A(1>»,35) 

A (15,1 •*> =QL/DQ 
F(.l>=Y(Ie>*-L< 

•*  ( 1 2 > = Y C -J  > 

*•  ( 1 3 > = Y ( 1 1 1) ) 

! (.U)=Y(Id)  + c i< 

O ( 1 5 ) = Y ( i lit 

no  ioo  ; = l . n x 

00  100  J = 1.4> 

100  F(i»J)=“  ( 1 ) * A ( 1 * J ) 

F ( * 1 . * « ) = M 

f (.4iL  II  s*<i 

* U ■ f 

1 i.J 
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FROM  COPY  FURNISHED  TO  DDC  


COMMENT  « 


COMMENT » 
COMMENT « 
COMMENT  I 
COMMENT  I 


700 

701 


702 

703 


70  4 


COM°:U*  FREE, DOUBLE  PM'.CISION 
SUBROUTINE  U0AT7  <X,NX,Y,N,PXY,°X0.7) 

SOLVl  FOR  ^nrK  + l)  ,PX0(K  + 1)  AND  X (K>1  > 

DIMENSION  a (NX ) , Y (N) »PXY  (N>  ,N)  ,PX  0 (NX, NX)  , O Fa  ( « • 5 I , OF Y ( 5 , 1 1 > 
JIM" NS  ION  Xl(5l,X2(f),cXYl(5»ll) 

OI  i-  NEICN  DO  (5  tll>  » DC1  (f  ,1 1)  ,DD2 (5,11)  ,DX  (5), 0X1(5), 0X?(5.) 
CA.L  JF7(JFX,CF'Y,X,Y,t  > , N ) 

»XY  A 40  3X0  AF  T SrNSITIVI7Y  MATRICES 
OJ=DUF'MY  MATPIX  = 0FF  I V A T I VE  W.RT  TIME  or  A SENSITIVITY  MATRIX 
USES  RUNOL-KUTT A METHOO  TO  SOLVE  OIFF,  i INS, 
f*fst  compute  srNsmv:TY  to  y 
OJ  700  u-l«NV 
JO  700  H = l.*4 
DO(L,M)=OFY (c,M) 

DU  700  J=1,N> 

D3(t,M)=DJ(L«M)+0FX(L,J)*FXY(J.M) 

DO  701  _ =1 , NX 
JO  701  m=1,n 

PXY1  (L,M>  =<-  J (L  ,M)  »T/3.  +PXY  (L  ,M) 

DO  7 02  u = l,N.< 

DO  702  M=1,N 
0.U  (c,M)  =JFY  (l  ,M) 

DO  702  J = 1,'4X 

Djl(L,M)=OOl(L.M)+rF> (L,U)*PXY1(J, 1) 

DC  703  l=1,N* 

DO  703  M=1,N 

3XY1  (u,M)  = (OJl  (L,M)  + DL:  (L,M)  ) ^T/6.  + PXY  (L,M) 

JO  7 0-  L S1 , NX 
Do  70a  M — 1 , N 
DD1 ( _, T )=JFY(l ,M> 

JO  70-  J=1,NX 

DJI (c,M) ajtl (l  ,M) +DFX (l, J)»  PXY1 ( J,M) 

DO  705  l=1,N, 

JO  70"  K=1,N 
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705 

706 

707 

70«* 

709 

710 

COMIt  N' 

500 

501 


•=X ; .,9>  = <JDlU.«M>*3.*PD<L»h>>*T/$,*PXY<L*MI 

JO  fuo  ..  = i » N a 
JO  706 

002 ( L * 9) =JFY (L  ,9) 

00  7 06  J =1 » NX 

ODa-(L.M)  = 0i  2 (L  *M)  + CFXCL  , J)  *PXY1(  J,  'll 
DO  707  u=l.NX 
Dj  707  9=1, N 

PXY1  (l  .9)  = t 30?  (L»9>*<*.-PP1  (L.H)*3.+0J(  L.M  > »*T/2,  *®XY  (l,M> 
00  70?  W-1.UX 
OG  7 0**  1=1,  \ 

3Jl<t»M=JFYU  ,M> 

00  7 0 a J=l,f> 

001(l.HI=0C1(L.MMDF>  (L.J)*°XY1(J,M> 

00  7 3)  w = 1 , *4  x 
00  709  1 = 1,  N 

PXYl  (L  »M)  = (00c  (L»  9 )**-.♦  001  (L  ,M>  +00  <L  , M )»  * T /6.  ♦ DX  Y ( L , M ) 

JC  710  l=1,NX 

JO  710  M=1,N 

F*Y  (L,M)=P>  Y1  (L,M) 

» CpI’UT1".  3‘  Nil '‘IVI‘rY  TO  XO 
DO  500  1=1 .NX 
DO  500  9=1, NX 
00 <u,  1>=0. 

00  500  J = 1 , NX 

J0(L,,'il  = J0(L,N)+0F>  (L.  J)*P>0  (J,M) 

501  u = 1 , NX 
00  501  1=1, NX 

3XY1 (L ,9) =00 (L ,9) MT/3. +PX0  (L ,1 ) 

DO  5 0 2 L =1 » NX 
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502 

503 


5 0“ 


505 


506 

5 Of 

506 

509 

510 

COMML  NT  l 

5 0 0 

410 

420 

430 

440 
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50<?  ,,  = lt  lX  mM  OOFY  FURNISHED  TO  liDC 

JD1 (w . M) =0. 

JO  5 02  J - 1 • MX 

DJI  ( _,  M = DC 1 (c  t M > *CF>  (l,  J)  »XY1  (J,  -I) 

J j 5 0 1 „ = 1 , NX 
jO  5 0 1 v.  =1 , IX 

pa  Y1  (l  .'ll  = <CD1  <L  *cr  (L.f>  ) *7/0.  ♦=*  0 ( L*r> 

JO  50“  =1,  N.- 

Dl  50“  P=1,NX 
J J1  <u,M  =0. 

DO  50“  0=1, NX 

1 (L  ,M>  «-TFX  (l  , J»*  PXY1<  J.M) 

Do  505  L = 1 * N a 

Jo  501  M^i, n> 

=xrl  U .mi  =(DC  1 (L»M»*-3,*DO(L*N>  > ♦ {/  3.  ♦PX0CL.M) 

J'J  5 05  0=1.  MX 
'JO  5 0o  M = l , n x 
D02  <<-,M> =0. 

DC  5 Or-  0 = 1,N, 

KD2  ( r ) = 052  (L  .M)  + CF>  (L.  J)  * ?XY1  ( J,  m) 

DC  507  L = t » N a 
0 J 5 0 / M = 1 . M> 

pxyi  (w ,m»  = or?  (l.m)*4,-dci  (L  ,m)*3.  ♦ojil.m n *r/?.  +pxo (l.m) 
JU  5 0.“  L=i»'JX 
do  50 j m=i,n> 

JJl  ( L » N‘. ) =0. 

Ou  50*  0=1, NX 

DJlfL.M|s3Jl(L.MI*tFXL.J)*P>Yl(J.MI 
DO  509  _=1.NX 
DJ  50  < M=1 , NX 

3aU  <_  ,M»  = ( JC2  CL  « M)»a  + r01  <L  ,M>  *00  (l,  I)  ) *T  /G,  *-P)  0 ( L . M ) 

JO  510  1=1,  MX 
JO  510  M = 1 » NX 
Pa0(L.^»=P>Y1(L.'1') 

IJPJAT!  ST  Alt'  V1CT(.K 
C A u L Gf>. A3^(X,f<X,Y»1  ,CX) 

DO  4 00  J = 1 » N> 

Xl<0>=CX<J)M/3.+X  (J) 

CAwL  G F A J 7 ( X 1 * MX ,Y , N » D > 1 ) 

JO  “10  0 = 1 » NX 

Xl(J)=(u\l(J»+DX(J))*'/b.+X(J) 

CALL  G r A D 7 (X'l.NX  ,Y  » K , C X 1) 

DO  i*20  J =1 , NX 

XI  <J>  = <DKl  ( J)*3,  ♦JX<JM«"/5.*X  ( J» 

CAlL  Of AJ7U1 , NX ,Y,f ,0X2) 

00  4 30  0 = 1 , Ma 

X 1 ( J ) = C DX  2 < J » * **  • -OX  l(J)*3.0X<jn*T/2,,x(JI 
CAwL  GFAJ7 (XI, NX ,Y ,N,0Xl) 

00  450  0=1, NX 

X ( J)  = ( CX2  ( J»  *4.  ox  1 ( J)  ♦OX  ( J)  >*T/c.  *X  ( J ) 

=_ruco 

END 
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ddofl*  3fftc:s:oN 

s jt  *ctr  :’i;  :»jv  ( a « 0 » k ) 

JlM-NilON  A (N, N)  ,n  <N,\> 
do  ii  : = i 
jo  in  j=i,'i 

10  3(.» j) =o . 

11  T(-.i) =1. 

;jct  a = i . 

no  ?o  k=i»‘. 

A<K=A(K,<) 

Of.:A  = J(  ■.  A*aK< 

JO  3 0 J= 1 1 N 

30  0 (K , J) =3(<»J)/A<K 

jo  *0  j=<,n 

40  A (<, J) =A <K, J) /AKK 

nc  5 0 k i = l . :j 

if  ( <.  i k:  ) r,o  ~c  r.  o 

OO  6 0 w = l ,(> 

t>0  a(<I,L)  = d(<l.L»-A(K:.K)*P(K.O) 

:f  ( k.  { iQ.m)  oc  to  'o 

K*>1=«>1 

DO  70  L=KP1.N 

70  A(<I,_)=A(M,L)-A(KI,K)*A(K,fc.» 

50  CO^TI.'#U: 

20  Co’O  JNU: 

X WKiTf  C6» 101)  r‘rj 

X101  FORMAT  (IX, "DC  TE^MlNANTs’M^l'-.n 

£ N J 
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COM'-;  s'  « 
COMsfUT « 
COHSF ST  « 
CO*  k ST  I 


COMVlnT  I 


CUMS:  ST  « 


‘•oil 


«•  TO 


*60 


i.  -<  0 

25? 
?$n 
?oo  i 


APPENDIX  D 

GENERATOR  SIMULATION  PROGRAM 

^iso.f  a nicHisf  s:k.'<.at:<:n 
jaij  :»  o data 

is::;a..  ftA*t  v<  ru ; cc*  c'U ' . ' p'jm  pa*- ams  t f *s 
ij  < ul.-.1  J US(  0 . Gf SL6i  CC-tY.'t  STb 

cot’:..,*  Fvtt.nouetL  cisicim 

COSMtS  •.''N.’-’oK  ,W 

3i  ir  No  K'S  1(4,  4>  ,Y  <1 1»  ./‘MO  ,?  (4»  , < <f*l  tAl  ( r’l  ,A2<$I  ,DX  (e  ) .3X1  (51  ,OX  2(5» 
Z Y ( 0 , 1 1 > , C 0 ( H)fC(4,lU  , j<  5) 
j-*-  a :/o  Fjiii 
0 s>I.  N F,,,jATASt" 

OP.  S " 

t-'.ar  in^u*  oa; a 

- £ F.  0 ( 0 » N,M  .t/.KFfT  ,W 
4 * 0(5  > . F‘  ,'1  1«<1»-1?»K2 

* l JO  ('■  I <Y(J),J31,M 
K 3 I 

**.  c X 3 • L 1 

COw  XSSMX.NX  ,Y  .N.Cn.CY.TM) 

* j=r/SiT 

ij  250  <3l,<r' 

lk3-:l2 
■* u k = ' l i 

:<3f*  <K-n 

JO  2 52  J<3UN2 

c-J-L  l>C*Y6  (X  ,N>  , Y , N.f/.CO.C.L'.  ?.  ’’A.TSJ 
7<3T<»* J 

ok,  \*F  A Jb(  X.UX,  V ,(■  , PH  > 

GO  *.5  0 J 3 1 « N' 

H 1 ( J > 3 i < < J > » ' J / .!  . ♦ X ( J > 

„al;.  o* Ajb ( hi , nh ,y ,s,r x l ) 

SO  4 fo  J 3 1 » N> 

Hl(J)s(PAl(J)*t'X(J)»»TJ/t.*X(J> 
uA.w  u» A Jt (XI .SX .V ,N,CX 1 ) 

DO  OO  J 3 1 . S/ 

Xl(J»3CMl(JI,(.*r'X(J>)*'J/d.*V(J) 

C A . L j*4i)rl  4 1 ,F(  *Y  ) 

Jo  •» c 0 J3i,s.- 

<1  ( J»  3 ( *2(  J>*  - .-rxl  (J)*.(,  ♦.'a  ( J)  ) *’J/2.*X  ( J) 

•OA.L  G.  a J6  (Hi  , NY  , Y ,s.px  n 
JO  uof)  J3l.N> 

a(J»3(Ia?(J»**..M'y1(JM2x(J>>*'J/c.*a(J» 
h ?(  J>  *.?•  ATS  (V  ( J I •>  1 ( J ) ) 

WO*  (10,2000  ( X 2 ( O . J3  1 . NX  ) 

: -m  jo -i  ( jx» 

XKiTC  (b./OOJJ  (7£'(Ll  » l 3 l »*  / > 

«"0,  (A*  (5<JU  ,t.) 

, SO 


Si!  ’ESI  warn  iucnuiu 

iwi  ix>]JOC 


a j;o°l . P>> . cn  1 1*' 

JV.^AY  2V.1 

iUC-OU'  ; X Xi'i'f  (X  . NX  , Y , «,C(tivY.'ll) 

COM 1-  MT  I C’i  40V  S*;u  Ot  i»IG  CO?  iM  F A*  AN  t Ti  > i.S:T:MATw  v 

COM-ILNII  0Y*~Af*.  0r>.  OF  >0  Ht  1 Y 

CljMM  IN  1 F < , - C N'  . W 

UiM\  ts  Ij'I  V (f.A  » . Y <f  ) .0Y  (,YX  ,N)  ,C1)  <Y> 

.1*1 
.2*2 
:s-s 
: •***• 

. ->*5 
1 o*  a 

: t»=H 
. ■«=•» 

*Y<;  i I*U< 

J*^-'  <2 . / J . > 

psj.  1-.1P  »2iv  S'  N.'a 
-’  = 2.  * J/J. 

DX=v  I.  *»>  • ( - *•  2+W*  • 2*  Y ( J 1»  * V (I  t)  I ) 

»(.1M,  FX»Yt:2)*r**2/^Y 

\ (.2)=. F\* Y I.LWr(:^»>Y(;.’)”2*Y(l()l*fFX‘  W**2/C1 

\ut»,K*Yi:n/Y(m*Ycii*Y<:ji*Y(.bi*fFK  »w**?/c>x 
\ <.•.»  = '*■  fk1  y <:  ?>  • r cc  > * «*w/r)x 
' (.F  ) s-_F<*Y  < ; ?»  * Y <T’>  » / O' 

x ini  . *1 . i 

Ou  101  J * 1 • N A 

ioi  oy(j.:i»o. 

.Y(;i..i*=*x{.n«‘tJ**2*Y(;t)*Y(T^)/.x 

or ( i l, ;2>*> c l >/y  < w » 

0Y(.t.:o)*-XC.  1)  *Y  <!  11  * Y**2*Y(im/0X 
Y(ll,.5l:^«(.llMl/'  - 
■>Y  Cl  1 . ! ■»)  = -X  r 1 ) /Y  (I>  ) 

0Y(;2..1)*w**.*Y<:?l*'2*Yt!t.)»*2*Y(I»»*tFK/TA*»2 
Y ( ! 2 ,.2»*-2*Y<:2»*v<:f>-’FN«w*2/jX 
)YC2.  :»»  3|  FN/Y<2-4» 

or  < :o. . i > *-  r < 1 2)  • • »•  ( fh»  w ♦ »/ox  *w*»i-*r«ii>,‘Y(;2)**2*Y(:6»»Y(im*i.rK./ox*» 
•)Y  ( 12,  . b)  ( :2  > *•  r*  Y Cf  > “h«  Y *1  F<»w*  *?/DX*»2 

3y i:  *-x  c *>  /y  c y) 

r ( : y.  : i)  * ■<•»-*  y ( r.’i  * y ( . y»  ■■  y ( if  t **  •* r ( m *i  cN/rx**2 
JY ( II. *2  I *-Y  C J> ' Y (IF  » * 2/OX 

ni;  u«-y<:2»»yi:*  )*.fk*w»»2/o> 

JY<; Y..-1  I'  F< / Y t . u ) 

JY  »:  1,.c)  s-V  C ’)  *Y  «:.!»•  t 2/0x 

• ♦H»»-*Y<il»*Y(;2»*Y(I.Y|*Y(A‘>)«Y(  . * 1 • _F  K/O  X * • 2 

jY<;.l,.3»-2*Y(;2)*Y(10)»Y(i6>»-*Yi:*J»*(FX*Y.**2/OX**2 

nci.nisM  r y > /y  <i“) 

i)Y(:*..iii-\(:ui*n*«?.Y(;n*Y(:mn> 

'yc-.:2»  *v  (.  - i/yc.’i 

oy(;-..:j»*-x(.<.i*(-i/y(;p»*w»»?»y(ti)*yu  ■»>  /cx  > 
jy  t:-.. .1)  *-x  j:  h • <- i/F  *2** • y (i <»woa > 

: y *«x  <:-»  /y  i;  »t 

j y » ; " , ; t » * - x i : * * * w * • ;?  • y y . t. » * y « r -r  > / ox 
jy ( . . 2)  ( .e » /y  t.  .?» 

' Y ( » ,%  • T 3 ) *-x  c;  r >•«*•?*  y ci  >*  Y U'Y*  / ’X 
or « r ir»  »*«:«•  »/y  c:7 1 

0Y(tv  , . A»  \ r.  01  •(-l/x  *2*'*Y  c Y» /0X> 

0 Y ( I a • I >>  =-X  Cr)/Y(!SI 
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COPY  JUKNlSKJiD  TV  LUC 

v.0Cti>  =-0*Y  Cfc  ) / (-  **2+W**2*Y  (Ifcl*  * 2) 

KtTU-.i« 

C N J 

COMMENT  TH  = INiTIAl  -CTJ?  £ N C-  L 1 C0  = G;-  * IVAT  iVc  OF  TH  WST  Y 

Th:4T:M.V»Y  CH/r) 

OU  169  M=1.N 
160  C0«M»=0. 

LO(ib) =***/(-  ”2+W»2'Y(>)*»2) 

COMPILE  u»CU?L<-.  D«fCl3;c*J 
j UP^jU  ~ i N i 0*  £Jb(X,N>  » Y ♦ N » u Y Y 
THIS  iUe-.OuriN-;  C l j LX/u’  = FCXCD.T) 

< (1>  =c£M3jA  P.  y (2)  =LAPr3p*.  F,  < ( J> = L £ 9RD A KO 

v U)  =•_  AM  306  9,  >(3)=lAPB0A  KO 

0O1NCN  :F<,r‘..Y,w 

Oa.^NSION  X (NX  >,Y(M  ,OX(M/>  , A(5,5>  ,F(5,5> 

12=2 

OAcL  AFb (Y,N,NX, A.F) 

00  100  1=1, NX 
Da (I ) =0. 

JC  100  J = 1 » NX 
3X(I)=CXC)*F  (I,  J)-X(J) 

OX  (121=  OX  C2)*fFK 
-i i J;N 
NO 

CC^:k:-.S  Fx  = :t''3U3LC  r?i  cision 

OV-'lAY  jVl2 

SU3r  JUT  IN.:  GCHYb (X,NX , V ,N,N2,C0 , C, 0, ZP,T, TH> 

CJ.1  tCN  iF<,^«K,W 

01  -IE MS  ION  \ (NX  ) » Y (N)  ,C0  (f.>  ,C  (NZ,N»  ,ZP(  NZ)  ,0(NZ  ,N>  ) , A (5,5)  ,F  (5, 5)  , 0,rH(4l 

• t * _ i c » _ o 
:i=i 
: 2=2 
: s=s 

14=4 

- 5 = 5 
.6=6 
i7=* 

; *=e 

J=3. 14 15  92o5  35  {.*79 

6*2, “P/3. 

S*-IV  T (2,/3.) 

CAlL  AFb(Y,N,NX.A,F) 

10  = 0 . 

JO  i»  0 J=  1 » 3 


COMNONI  I 

c 

c 


100 


125 
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CO 

10  = 

I 3*  A < 

; l 

, J»* 

*<  J) 

*Q= 

0. 

0 0 

50  J = 

•» , 

c 

50 

i a* 

: }♦*< 

i •» 

. J)  * 

X(  J) 

GO 

•io  : = 

l. 

N 

00 

40  J = 

l. 

HZ 

90 

L ( J 

=o 

• 

00 

95  1 = 

l. 

N X 

00 

95  J= 

i. 

hZ 

95 

D(  J 

. i)  =0 

• 

DO 

i no  r 

= i 

, 3 

0 <i 

1 .*  > = 

0* 

A ( * 1 

.1)* 

CO 

i ( w*  * 

* T H ) 

D(. 

2,.  ) = 

3* 

A ( I 1 

,1)  *• 

COS  (WT  + 'H-3) 

0<- 

3,1)  = 

o* 

M (i  1 

, 1 ) - 

CO 

SAW*" 

* 7H*3> 

100 

o(ic,:»= 

A ( 

*2,1 

) 

:o 

105  : 

= C 

.5 

0(1 

i,n  = 

‘AC 

c , * ) 

T ♦TH ) 

0 < 1 2 » . ) = 

-S 

»A(i 

H,  I ) 

IN  <W* 

T *TH“ F) 

105 

.)<: 

3,:»  = 

_ ( 

'AC 

-,i) 

♦ c 

I 'i  ( M* 

7 ♦ TM*  =) 

0 0 

nn  : 

= 1 

. u 

Z?< 

: > = o. 

00 

1 10  J 

= 1 

. 4* 

110  P(ll=Z5(I)*i)(!*J)‘X(JI 

00=  Y <:i>MV(:-»*V(IJ)-Y(15l**2l-V<I2)MY(IJI*Y(I2l-Y(I5»*Y{IJH 

• •vciJ)Mr(>i*Yi;i)-rcz)*y(:5)i 

CQ=’Micl*YC7)*Y(;7)‘*?  * 

01  = A(:i.;i)*^C:iM-A(ilti?>*X(I2)+A(U.I3)*>(I3) 

oci.in=-A(.i,ii»*oi 

v.(.l.:2)  = (-Y(I  3)  *>’  ( I 2)  *Y  C5)*X  (13)  1/  DD*2* A (I1»I2)*D1 

c«.  i «:3»  = m;<.  >*x  < 1 1 » - v c2i*x  < 1 2 » *r  <:s  >*x  <:2>-y<i*o*x<i3>  »/oo 

* -<-(:if:i»*2*4-uitH>)*oi  ‘ 

C(*l.:‘.l  = (YCX)*V(ii)-Y(i?)*x(I3)l/DJ-A(I2,I2)*01  t 

C<:i.;5)=<-2*Y(*5)«X<ll)*Y<I3)*X(X2)»Y(:2l*X(I3ll  /OG-2*  d <12,131 *01  . j 

GC.l  ,lb>  =-A 

cc.i, :')  = (•  <;o-x  <;?>)/oq*  <4<*5t:*3>«-2*Mi*.:5>  >*oc 

•00  1 20  *=1.5  j 

ci*3.*i  = o*c  ci»i)*rxb<w  ; 

t <i’.;)*i*Li:i.:)*cos(w*T-fTH-=)  I 

120  ci;i,.i=>*l(;i.:)*co£<h»t+-h)  : 

i io  ; =u. 7 

o <: i.:  i =-i*cc  l.  :»*sif. <h*r  + TH*=»  I 

cc2..»s-i»on  l. ;>* ?:►. (w* t^th-c) 


130  - (.1 >=-$' C C 1. I>*  olN <9«t*tM) 


o o 


. ■ 
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ci.mii  = (»ini,)((i?i-r(:ei*>  C4* i/co-aii i,m,02 
c(;w:?)  = i*r(.M)*v(ii)4rini*xiiji  )/oo-2* a <;i. t?> *02 
n:i.,;ii5iir(’.i;i-v(i2ii*s(ii)KV(iii*?,if(iJii,*(i?i*vi:?i,)(iiJii/op 
-cc3,;3>  *2*;  ci.:  n>*32 

=--(i2,I2)*32 

cu,:r,)  = (Y(:3)*>(:i)-Y(:i>*x(i3))/o?-2»A(i2,i^)*u? 

DO  30  1=1,  tZ 
JJH(l>  =0. 

oth ( :n  = -s*  ( : *s:n<w*t ♦thi  *:q*cos («i*t  + tkh 
3r-c2)  = -r*(:c»s:N(w»T*TH-^)*iQ*cJs<H»r*TH--n 
0TH(A3)  = -i*  CC  *SINM*’*TH+C)  ♦IQ*C0S(rf*  TtTHt3|| 

JO  33  1=1,3 

C (-. Ib) =0  (L  , It  >*0!> (L) * CO  1161 

c c 2 1 = c c . if ) *crH  a j *co  c a > 

*.r  : J U 
- NO 
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i<  J-UdLr  er»ClSIOW 
SOd5..JTlNl  AFt <Y,N,NX,A,F) 

corntn.'i  *h.s  sis' dd”:*.;  cc^-’ot.s  a=:n»mu  and  f=-.-<»-a 
Coitjn  tFK,n.K,W 

o_-  \jiwN  y cn  .a  (n>  ,,u  > , mnx  ,nx  ) , ms> 

* 1=1 
.2  = 2 
. 3 = 3 


I 6 = n 
.7  = T ^ 

. d = » 

* *=** 

- u = in 
.11=11 

DO  3 0 1 = 1, -.X 
00  3 0 J=1,NX 
30  A(J,:)=0. 

00  A n . = 1,.\'A 
DO  -*0  O = 1 , N X 

WO  F ( J . 1 ) =0  . 

01  = Y (.  -,)  *Y  (13)  -Y  ( I 5 ) * • 2 

J 2 = Y (13)  *Y ( I 2 ) - Y ( I 5 I * Y (13) 

0 3=yC‘-)*y(;3)-y<:?>*y(:[t) 
u=Y  < i 1 ) * Jl-Y <: 2>  ♦ 02- Y ( 1 3>*  0 3 
D ) = Y f. t>)  *Y (.7) -Y (17) *• 2 

a (.1 n =ji/; 

A (.1,12) = - l 2 / 1 
A (. 1, ’ 3) =-03/G 
A C 2,  I 1)  =A  ( 1 1,  12) 

A (. 2ii?) = (Y (.1 )*Y ( 1 3 ) * Y ( J 3 ) * * 21/0 

A(.^.:3)=-(Y(11)*Y(I5)-VC2)*Y(I3I>/J 

A(i3.ii)=A<:i,;3) 

-»  ( . 3«. ?)=«(. 2,13) 

A(.5,.J>=<Y(.1>*Y<1-I-Y<.2)**2)/J 
A ( 1 w . 1 ■»  I =Y(17J  / OQ 

> /oo 

- d'. =a  ) 

A (I;  . . r I =Y  Co)  / JO 

- ( _ 1 > = v ( i ,>)■*■-  L < 

( . 2 1 = y « : j ) 

- ( . 3 ) = v ( : i o » 

- (1  * > = Y ( Zn > ♦- 1 < 

( . 5 ) = Y { . 11) 

,>J  too  1=1.  NX 

o 100  0=1  , 4 x 

100  c (.,  J)  =-' (.) * A (I.  J) 

F(ii,:-i =w 
-*<. •*,:  i)  =-w 

1 t T >J  %.  N 
' 41 
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COMiLNTl  6RAF1.FP 

C t-  r.  A D If.  JP  T D 5 CC  LUf-AL  OF  (.ill  DATA  A NO  Pw3T  THP  CN  CIJIVIlW 

C HI'  ANY  <t  Y K>  SLl  N-YT  PLOTS 

C inP  = Nw.  J F POIN’3  ctP  G-APh,  i43  = NO.  vF  GRAPHS  Cl)c  TO  10) 

C S=A9  AY  JF  SCAll  F A CT  C c 3 FO  -ACm  GRAPH  - 3 < J > =L  4)  Gt  3T  NJ.  ON 

c 

Oil!  NS  UN  M5.-5l).:>C5.4.51),iY<5.451>  o ( 5 ) . J>  U51 ) 
jor  4 5 , "O  A'  'i*.  A , Lf>  =999  » AT  T =**I  ** 

•V£AD  F*-t.  C5>  f ",  NG 


At  AO 

F f L 

t 1 

C5 

) 

CSC J) , J = 1 .NO) 

1 

CO  1 

9 < = 

lc 

») 

p 

"l-0 

F<E 

i 1 

15 

, 

• f 3=900)  CxC  J.K) ,J=1.  43) 

DO  1 

) J = 

lc 

> N 

G 

19 

I A ( J 

,<)  = 

A 1 

CO 

, 

<) * 100/S  C J) 

J=1 

20 

CO  NT 

IN  Cl 

C All  '•■SCO) 

C Aw w SHOO’. 

ca.w  *c_o;co,o,o) 

CAwl  TFLTT  CO, 250. 2) 

CAw-  ’PLOriSOO, 250,1) 

CA_l  T AH. SCSO, 125, 45 0.200, 10, 10, 0,-1001 
DC  JG  <=1.M» 

30  jx  (<)  = :<  (j,o 

CAuL  TPOIN(Ja,1,,I«F, 1,1,0) 

J=J*1 

IFCJ-NG)  50,50,60 

50  CA^.  VA*IS(50. 375. L50, 200, 10. 10,0,-100) 

CO  4-0  K = 1.NP 

40  jack>=:<  c j,o 

call  TPOiNCJX, l.,NP,l, 1.0) 

J = J*1 

lc<  J-f.G)  70,70, 60 
70  C 4_l  1M3JK0) 

PAU3  * 

60  r I 20 

60  C Awi.  Af'JOi  ( 0 ) 

PAl;3‘ 

G3  T i;  1 
900  STjD 

9A9  3T„  = “ :-'CF  IN  0*1  M DA  76'  J" 


-OT  j 


n r.  n r.  o r.  r. 
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APPl.NDIX  E raUM  COPY  FUKNlyHKL)  TOUDC 


STANDARD  PARAMETER  ESTIMATION  PROGRAM 


>•  pMcru  is  PAH»Mf  rtH  towoH  fac  row  of  initial  uuess 

>*  lif.LFAC  IS  FOR  I HP  I NC  HE  ME  N T OF  PFACTR 

► ♦ NINPAM  The  NIIMIILR  OF  REPEAT  OF  PF  AC.  I R *1>ELF  AC 

► ♦ NRCPIT  IS  NOMrtf P OF  ITERATION  OF  PARAM  ESTIMATION 
E*  NSAMPL  IS  FOR  THE  NUMBER  OF  samples 

► * NIRUN1.E  IS  The  number  of  rungf.  routine  repeat 
PE  AL  T R = 2 . 0 

oe i r ac=- o. a 

N I NP A R - r 

nrep I t = I a 

N IRNGE-  I 
pi  imi r * c . 2 s 
ML  -0  . 

UME GA= I . 

DEL  T- 0. 2 
N V-S 

DIMENSION  l l MT  X I •>  I . MMMT  X (5  ) 

0 I MEN  SI  UN  SSA l(S) 

DIMENSION  XR(S>.AR<S.S>.U(S>.VR(3>.CH<3.S) 

0 1 MEN  SI  ON  X(SI.A(S,S).Y<3).C(3.S) 

D I MEN  SI  UN  S X I S • I I ) »AP(5*  1 1 ) iSTI.li  I I ) • C P < i.l  I ) 

I)  I MEN  SI  ON  XSI  S)  ,US(  S > . W I NV  ( 3 .3  ) ,ST  T ( I I ,3  ) 

DIMENSION  Svr«(ll.5).SUM<ll  ,1  I I .ISUMU  | ,||) 

DIMENSION  VD I O.VWII  I I.TYTKII  )• TSUM I Nl  I I . I I ) 

DIMEnSIDN  Df  LP(  I 1 ) iPNl  I 1 ) 

DIMENSION  PTMUI  ( I I ) .PR A T I Ol  I I ) 

DIMENSION  LMTKIIl).MMTX(ll) 

DIMENSION  ESTMATIIS.il) 

COMMON  AR.A  » U , O S 
u<  I ) -0  . 
u(  2)  =0. 

ui  n = sum  i < )op  i 

U ( A ) - 3 , 

OI  •>)  -0. 

Ml  ADI  ).  10)  PI  .02  . P 1 .PA  ,PS ,P6  .P/ .PH ,PR .PI O ,P I 1 
I l)RMA  I ( I F l 0 • S / Jf  | 0 . •»  ) 

MM  I H (h.lOA)  P 1 . P,’  . P 1 .PA  , PS  ,P6  , PT  . MB  ,E>  3 . P I J , P I 1 
E CIMMA  I ( l HI  , l 0 X . 1 1 Ml  » PAR  A ME  TER  V A L UC  S • ✓ / / I I 0 X , 1 1 | 0 . S ) ) 
P I RUE  { I ) -P  l 
P I M ol.  I 2 ) =P  2 
P I MUl  < I)  I 
•IMOE"  ( A ) -PA 
P I RUI  ( S I -P  '» 

P T MUl  ( F I -PE> 

P I RUE'  I / I -P  / 

I'  TRUE  Iril  -Prt 

P T Rut  I ) ) -P'» 

PIP  UE  III')  =P  l 0 
P TRUE  III)  -P  I I 
IMI  30  U I 1 1 . M 
DU  30  0 3 - I , S 

A I I . J ) s i . 

AR  I I . J ) • 0. 

DU  II  0 1=1.3 

DO  110  1 = 1 .S 

Cl  I • J ) = 0 . 

CR ( I . J 1 =0  = 

CO  = -P  1/P  I 

COKL3  = P 1*  (P2-P  | ) / (P  1 * ( P2-PA  ) ) 

CDF  =P  1*  ( P I - PA  I * I P J-P2  ) / I P I ♦ I P.1-  PA  ) * I P2  • PA  ) ) 

i f *CO  ♦C.DKO*  I I . -PA*  (P  J-P4  ) / ( P J*  (P2  - PA  ) ) ) 

l E O = ( P I - PA ) • < P 3-P* ) / ( PI • I P2-PA ) ) 

f.F  KO  = (P.‘  -P  I ) • I P3-P4  ) *PA/  ( P I * ( P2  -PA  ) • *2  I 

LX  013  = ( P -PA  ) /P  I 

( X O =-  E*.»  /P  I 

(XDI  -IP  l-P2)*P4/IPI  • IP  1-  PA  ) ) 

LU--R)|/IO 

iUKUl  S()M  Mi.  O*  *2*  CO) 


U I MINS!  ON 

0 I MEN  SI  ON 
D I ME  NSI  UN 
D I MI  N SI  UN 
D |Mt  MSIi,N 
DIMENSION 


PKSCXQLNO  FAflS  bLAMC 
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A(  I . I > - < |P|0*BU/PM*CO 
A(  l-JIOHL)/(>M*CD*U 

A ( I , J » - < l>*|0*HL)/PI|«CDf 
A ( I . A | -UHL  .A 
A<  2.  | >-(  l./PStACKni) 

A ( 2 • 2 ) = I I . /PS  » »CKO 

a < 2.  1 t = i i . /ps » ackof 

A ( ) . 1 ( --  ( I . /PS  ) #CF  O 

A ( *.  2 I = ( I . /P6  ) ACF  KO 

Al  1.  t ) = ( I . /PS  ) ACF 
A < 4 . It-  -<)Ht'UA 
A ( 4 » 4 ) -I  ( P I 0 ♦ RL  t/PHIACU 
A ( 4 , jl  - ((PI  OAWl  I /P8  I ACUKU 
Al  S.A  t = ( l./Pt)«f  QKO 
A ( 5 .S  ) = < I . /PI  » »r.Q 
C ( I , l ) - ( | . /P  ))»CU 
C(  l./IMI./PJ  )«COM) 

1(1.  I)  -<  I . /P  I ) • CDF 
C ( 2.  I ) - ( l . / ( P I I *P6  I I A ( -CL  O ) 

C ( 2.2  ) = ( I . / (P  I I AP6  ) 1 A ( -l  FHI)  ) 

c < 2.  I ) = < t. /(PI  I *Pft I I a (-Cf I 
C(  1.4t=(l.  /PH  > AC  0 
C ( 1 . S ) = ( I . /P«  I AC  (in  U 
DO  ."Pj  ( = 1 . •» 

III)  2 IS  .1=1  .s 

2A>5  AW(  I.  J)  =A  ( I . J) 

OU  2S»tJ  I =1  . i 
OO  2 R6  J=  I « s 

2<#<>  cmi.ji-cii.ji 

i*R  I Tfc  (6.  101  )(  (API  I »J)  . J = I .S)  »l*l »S) 

101  FORMA M I HO  . 1 OX. • NT  X A R • / / / ( I 0 * . SF 2 O . Ti  > 

»W  I TE  ( 6 . 10.*  > (ICR(I.J).J  = I.S).I*».JI 
1C 2 FORMA T(  I HO.  I OX .* NT  X CR,///(l0X,sF20.M» 

C 

CCA  AA  ♦ AA  AA  A AA  A AA  A I NI’Ul  NO  I St  IS  ADOE  O A A A A A A A S I AH  T 

c 

CALI  M1nV(A.S.OD,LLM1X.NMMIX> 

C AL  L uMPPl)  < A . U.  > SA  I . S . S . I > 

no  i - i • -> 

2M  t SSA  I ( I ) - -SSA  I ( I ) 

«W  I It  (6. 28 2 M SSA  I ( I t .1  - » . I 
20/  FURMA  I ( 1H0.SF  20.  S) 

L)  i NPU  1=0.01 
( UMAU I =SSA  I ( I I AO  I NPU I 
COMAU  2=  • SA l( 2 > AO  I NPU  f 
CuMAU  ( = SSA  I ( J 1 A 1 | NPU  I 
CuMAUA  ’-.iSAI  ( 4 ) AO  I NPU  I 
C(.MAUS  = S SA  I ( S ) AO  J NPU  I 
CGMAU  l-AHS(CUMAUl  ) 

COMAU2  - A US  ( f.  L.MAU2  ) 

COMAU  )-ABb(C(»MAlJJ| 

C(.MAU4  - APS  ( C0MAU4  ) 

CCjMAUS  = ArtS(  O (i  M A U S > 

1)0  2 OS  I = I . t 

00  2 OS  J - I . » 

2 0 *.»  W tNV(  I . > > = 0 . 

A I N V ( I,  11-0.02204/ 

4 1 N V ( / ) ••  1.  H41) 

m IN/(  i.  1 ) =0.02204  r 
A I NV<  1.1  ) - 1 . 0 
«INV(2»2)=I. 

« I N V ( 1.  1 ) = l . 0 
III)  PPrt  NP  ah  a M i | , N 1 NP  am 
•>'l ( 1 > IP  I out  ( 1 ) • P F Al  1 R 
•N  ( 2 I -P  I RUI  ( 2 ) APF  AC  I R 
• >N  ( » ) -P  ( RUt  ( 1 > API  AC  T >1 
■ »N  ( 4 t -P  I RUt  ( 4 ) APT  AC  I R 
»•.(',)  =>»  r PUT  ( s ) APF  AC  1 H 
»*’.  I h ) Ip  I RUt  ( t)  ) AI»F  AC  TR 
.'N  ( ' ' P I PI  IF.  ( / ) API  Al. 

'■-.(Ml  -P  1 UJt  ( H t API  Ai  I R 
PO  III  -PI  RU(  ( DAl’l  »(  ’ » 
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i>ni  i o i -i»  irue  t i o > • pf  ac  rn 

PN(  11  I =P  (RUE!  1 I > *PE  AC  tH 
IHN04RAM.I.T,  7)  r, Q TO  II 
PN  I I I =1*  TRUE  I I ) • I .6 

p.* « 2 ) =p  truf  < a i • i .« 

PNt  |)  JP  f HUE  ( J ) * I . 2 
PNt  4 ) =P  I HUE  ( 4 ) • I . fl 
t'N  ( | :RIHUfc'(S)*O.S 

>*N  ( 6 ) = P TRUE  (61*0.7 
PNI  7 1 -PI  HUE  ( 7 1*0.6 
pn(  a i =p  true  < si*o  . n 

PN( 01 =P I HUE ( A)*0.4 
An  ( I 0 I TRUE  ( l 0 ) *2 . 0 
Put  1 1 ) -P  IRUE  ( I l ) *0. 4 
i i on  j (-1.11 

PHA  r IO(  11 =PN I I > /P I HUE  ( I ) 
I CONTI  NUl 

DU  4 1=1.11 

4 ESTMA  T(  I .1  l-PRATIOl I ) 

DO  9<i9  N I T a I • NRE  P I T 
t * t =5 
U2;S“) 

( A 3 =‘>‘3*1 

l X 4 =5  t 56 
I *5=66555 

I ON  =9  V I 
IFN  = 3J  3 


1 ON 

= 777 

DU 

20  l 

1 = 

% 

s 

<P  ( 

1 » -0 

• 

X<  1 

1=0. 

DO  2 02 

i 

- 1 

5 

2.  - 2 

A ( 

1 . J ) 

- 

0. 

DO 

20  » 

i 

= 1 

1 

1 

5 X 

( 1 . J 

) 

-0 

2 C 1 

AP  < 

1 . J 1 

- 

0. 

oo 

20  1 

t = 

. 

.1 

DO 

20  4 

t 

= 1 

6 

204 

( ( 1 

. J > - 

c 

• 

OO 

20  J 

) 

= 1 

1 

1 

f.P  ( 

1 . J ) 

- 

0. 

2 0 i 

5*1 

1 . } > 

- 

0. 

Dll 

206 

1 

= 1 

1 

1 

DU 

2 06 

J- 

. 

1 

T ■; 

DM  ( 1 

• 

) ) 

- 

0 

• 

2 J 6 

»OM<  | , ) 

) 

= 0 

ui ; 

209 

I 

= 1 

1 

1 

2 C 9 

1 Y A 

(11  = 

0 

• 

P 1 - 

PN  ( 1 

) 

P.»  SAN  ( 2 I 
P 1 sAN « J | 

P 4 -PN ( 4 | 

AS  = PN  ( 6 I 
>•6  = PN  ( 6 I 
p r san < r ) 

I'M  -HN  ( H I 
P * -PN  < 9 I 
P I C -Pul  I 0 1 
' I I = >»N(  I I > 

>.U  s-P  J/P  | 

'.DMJsAJ*  IA2-P  | ) / IP  I * ( P2-P4  1 I 

CUE  -P  1*  ( PI  -P4  > • ( P3  - P2  I / ( PI  * ( P I -P4  } ♦ ( P2  -P4  I | 
! I =00  ♦(.UKi»*  ( l . -P4*(P3-P4  )/|PJ*(P2-P4  ) ) I 
If  D=(|i|-A4  ) * (P  3-P4  ) / (P|  * ( P2-P4  ) I 
V TKUMPf  -PI  )*(PI-P4)*P4/(P1  * (P2 -PA  I **2  I 
(.  » OO  s |A,’-A4  ) /P  l 
(.  K 0 = -P2  /P  l 

C A PE  ( P J-P  2 I *P4  / I p | * ( p 4-HA  ) I 

l. vi  = - P fl/P  1 
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( OK  U i 'id*-!  1 < C o*  • «*  ♦ c 0 > 

■\  ( I . I > - < <P  l Of  Pi  1/PJtKO 
4 1 I . ) - 1 ( P 1 0 ♦ ML  1 /P  J l ♦COKI) 

At  I . i)it(('IO*»l.)/PJ)*COt 
A ( 1.4)  OVI  oA 
X ( .*  . I ) ; 1 l . / P S 1 A C KOI) 

U^w’l  l l.'filKKl) 

A ( 2 . J > - < I . /P'S  ) AC  KOI 
A ( 1.1  ) - ( I i /Phi  *OrD 
X < 1 . i ) = ( I . /Mfc  1 *' F K ) 

A I 1.  1 ) - « > . /Pf>  ) *C  F 

A ( A . I)  - -l)  ML  u A 

a(  a.  a j ( (pi  o»m  ) /*j«  i •<  o 
A ( A , 'll  - ( ( P I 0*  ML  ) /!•«  ) *'-0KU 
A(  «=>.  A » = ( I . /PM  I *C  JKO 
At  6 .>  > = < l./P1)»C(l 
((1,1)  II  ,/Pl)«Cl) 

C ( I . 2 I - < I . /P  1 ) ACOKO 
CM,  11  - I I ,/Pll»COF 
l.  ( * . I 1 = ( 1 .✓«-»»  I AP6  I 1 • ( -Cf  U 1 

C ( 2 . «?  ) - < I . / (PI  I AP6  | ) A ( -CF  Kl)  ) 

( 2.  1 I - < l,/(P||Al>6l)*(  -CF  ) 
f ( • . A ) - ( I . /P«> AC  Q 
C ( •.  i ) - ( I . /P8  ) AC  OKU 
t.OF>  I sPI/i’IAA? 


KwST  QUAL1TY  PRAcncian 

K i)  UhM;iHJ£D  TO  UDC 


CLP2=C. 

COM  U-t./Pl 
C DP A = 0 « 

COKOP  1 A -M2 AM S/ 1M 1A  A 2 A (P2-PA  > I 
COKOP  2*M  JA  ( P l -PA  > ✓ ( P 1 A (P2-P4 ) A A2 ) 

COKIjP  1*(PJ-P|  I / < (P2-PA  ) API  ) 

COKOP A = < P2-P 1 > AP  3/ JP  » A (P2-PA 1 A A 2 I 

COFP I =P  JAM AA  <P  J-P2 I / ( < P J-MA  > A (M2-P4 ) AMI  A A2 I 

CCJFP?  -M  JA  (Ml  -MA  > / (Ml  A ( M2-M4 ) A A 2 > 

C l)F  P J *(  ►’  l - PA  ) A (P  JA  (P  J-PA  )*PAA(P2-P1)I/(MI  A(  M2 -PA  |A(P  1-MA  I A A 2 I 
Cl1!  I’A  ' P IMP  1-Pi  ) A ( (PI  - PI  ) A ( IV-t’4  I * ( P 1 - P?  ) A ( P 1 — PA  1 | / 

• ( P I A ( (M  I -M  A ) A ( M2  -1>  A ) ) A A 2 ) 

CKOOMI  (l>A-P2)/P»AA2 
C KOOM  2 - I , / M I 
!"K  POP  1 = 0, 


(.  K OOP  A -=  v I . /Ml 
1 KOM  I =P  /M  ; A A 2 
i klUV-l./Ml 
OK  OM  J =0  , 

C K L)P  A 0 . 

O'  orp  | = ( M2-M  U AM  A / I ( P 1-PA  > API  A AP  1 
CKOI  P 2 =M  A / ( (PA  - M 1 ) AP  l ) 

CkOFP  J = i»  A A (P2-PA  | / (P  I A (P  3- PA  ) A A 2 ) 

CK  Of  P4=MJA(M  1-P?  > / (M 1 A (P J-P4 I A A2 I 
CfOMI  = ( P 1-PA  I APA  / ( (P2-PA  I API  A A?  1 
t F OM2  - (M  l - PA  I A (MA  - p.i  I 2 (PI  A ( P?  -PA  J A A 2 ) 

C f OM  J = ( Ml -PA  I /( PI  A ( M?- P4  > ) 

CF  UP  A < M A A | M2-MA  1 A M2  A (PA-P3  ) API  A ( M 1-H2  I I / (M  1 A(  P2  - MA  I A A 2 I 
CFK0P|3(PA-P3IAP2APA/(P|A(M2-PA|AA^AM1 | 

IF  Kl»‘  ,»-•»  A A ( P I -PA  ) A ( 2 . AMI  - M2  -PA  > / ( P l A ( M2  -M4  I A A J 1 
CF  KOM  1= ( P2-P I I AM  A / (PI  • ( P 2 - P4  I A A 2 1 

Cl  KOP  A = < P2-P I | A ( P2 A ( M 1-PA  I A PA  A < M I - IV  ) > / I M I A ( IV  - MA  ) A A 1 ) 
CFP  I -ail*  I ♦ COKOP  | All.  MA  A ( M 1-PA  1 ✓ ( P 1 A M2  — M JAMA  1 ) 
CFM2=LDM2Ft  DKDP2 A ( I . -MA A (PJ-PA  > / (P  1 AP2-P 1 AM4  ) > 

• ACOKOAM  AA  ( P 1-PA  > / (P  1A(P2-PA|AA.») 

CFP  1 — C OF*  !♦  COKOP  1A  ( I.  -M4A (P  J-PA  > / (PJ AM2-P  1 AMA  I I 

• -f OKU am A A A ? / ( ( M ? - PA  I AP  J A A?  I 

CF  PASCOMAACOKOPA  A < I . -PA A ( M J -PA  1 / (MJ AP2 -M  J APA  ) ) 

• -CDKO*(  (P  1-M4  l*M?F( P4-P2 ) APA ) / ( P J A ( P2 -PA  ) A A2  » 
t OP  7=M8/P /AA2 

C0PH-- I . /P  f 

C 0 K QP  7 » ■?  , *j  A C Q A ( 2 • A PB  / P 7 • A 2 - I » / P 7 l / C OK  O 
O)KQP«s0 .5 A ( 2.  APfl/P  FA  A2-  | . /p  7 » /COKO 
10. 

(U»  10  I NKK  ~ I , MOO 

t-O  60  C - l . nIMfaOF 
OO  AO  J-I.ll 
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THIS  PAGE  IS  BEST  QUALITY  FRACTICABU 
FHUM  001  Y i'UKNlSHJED  IV  HOC  


0(1  4 5 (<1.5 

Ub(  I | >At'l  l«J) 

AS  MSI  i > *S * I I . J) 

CALL  MUNGE  ( I.Oa  J.AS.NV,  *> 

00  SO  I * I .5 
SO  SXI  I . J I = XSI  I I 

AC  CONTINUE 

CALL  GA  US  S I I M >Cli  MAUI  . 0 • • V V l t 
CALL  GAUSSI I X2.CGMAU2 .0. . VV?) 

CALL  GAUSSI I X3.CGNAU3.0. .VV3I 
CALL  GAUSSI I X4 .CGMAUA .0. . VV4 I 
CALL  GAUSS! I XS.CGNAUS .0. .VVSI 
U!  I I =VV  | 

u<2i>vv; 

Ul  3) * V V 3 

u(  a i =vva  - sum r < i.  i 
U(  5) =vvs 

CALL  MUNGfc  ( T .OEL  T , XH ,NV  . I I 
U< I ) -0* 

Ul  2 I *0. 

U<  Jl-C. 

U!  A I * - SUM  T ( J.  I 
Ul  S > =0. 

CALL  HUNGE ! T .DEL  I . X .NV .2 > 

AP  I I . 1 I =ICDP1 A x|  l ) KOKOP  I * X I 2 I ACDFPI AX  I 3 I I A I P I 0 a ML  I /P 3 
API  1.21  = ICOKOP2A  XI2  > ♦CUFP2AXI J|  ) A | P | O ♦ RL  > /P 3 
API  I . J)  - I COP  3A  XI  I ) ACOKDP  TAX  (?  ) KOf  PJAX  I3))A|P|0»RL)/P3- 
ICDAXI  I I ACOKOAXI2)  AC  l)F  A X ( II  I A | P I 0 A PL  I / P I A A2 
API  I . A ) .-  ICOXOPAA  X|  2 ) ACOFPAAX  | I ) » A ( P I 0 ♦ Ml  I /P  3 
AP  II  . I 0 > MOOA  XI  l 1 ACDKOAX  I?  > A(.Of  AX  I J > >/f»3 
AH  I 2 « I I - I C KOOP  I A X I I IACXOPI  A X {/  I ACKOf  PI  AXl  ])  )/l»5 
AP  I 2 « 2 I (CKOOP2AXII  ) A<  KO m2  A « (?)  ACKOA  P2  AX  I I)  >/P5 
AP(2,  1)  CKi)FPJAX(  ll/HS 

API  2.  A > - ICKOOPAA  X I I I AC  KUA  PA  A X I J > I /PS 
API  2,  --  I CKDD  A X I 1 » AC  XI  > A X (2  I Af.KOT  AXl  < ) ) / PS  A A 2 

AP  I J.  I ) = 1C  FOP  I A X I I ) ACf  KOPI  A X 12  ) AC.EPI  AX  I I I I/P6 

API  3.  2)  - ICFOP2AX  | 1 ) ACFKOP2  AX  | ? ) aCFP2*X  | I)  ) /PIS 
API  ).  J)  t ICFlJP  IAX  I 1 ) ACA  KOP  TAX  I 2 ) ACFP  IAX  I 3 > )/P6 

API  3.  A I - ICF  UP  A A XII  ) A CF  *DP4  A X | 2 | ACT  PA  AX  | 3 ) )/  Pft 

API  I.G)--<CFOAXtl)ACFKUAX|2)ACT  AXl  IM/P6AA2 
APIA.  / ) - I COP  FA  X<  A | A I OK  OP?  A X I -»  ) ) A ( P I (j  * ML  I / P« 
A"|A,«)-|PI0aML|aICUKQPH<'PH-lCKU/PHA*2)AX|SI 
API  A*  10)  -ICOA  XIA  ) ACOKOA  i (SI  l/PH 
APIS,  t I - If  OKQP  7A  A ( 4 ) AF  UP  / A X I *>  I ) /PU 
AP ( 5,  rt I I 0 OKUPrt*  X ( 4 ) AC  CPU  A X ( S ) I/P  I 
AP(i,U|--ILUKQAX(4|aC0AX(S)>/P#AA2 

1 I • f)  t l * 

60  CONMnUA 

CPI  l . I I -AP I I . I I / (PI OAHL I 
CP(  1,2)-AP(1.?)/(PI0a  ML  I 
C P I I . J I - AP  I I . I I / I P I 0 A ML  I 
CPI  I . A I - API  1 ,4)/<P|0AHL> 

CPI  2.  I > --API  I.  I I /PI  I 
CPI  2. 2)  --API  1.2) /PI  I 
CPI  2.  J I - - API  .1,  I)  /P  I I 
CPI  2.  A ) = -AP I l.Al/PII 
CPI  2.  ♦>  > --API  1.6)  /PI  I 
< P I 2 . II  > -C  P I 2 . 6 ) A | P6  /P  I I > 

CPI  J.  n APIA,  n/IPIOAKU 
CPI  J.H>  - API A.rt)/(P10ARL) 
oAlL  uMPRO 1C . SX. SV . 3 .5. I I I 
on  65  l - 1 . J 
UU  AS  J * I . I I 

AS  JY I I , J > -SV I I , J I ALP  I I . J) 

C Ac  l >jMi’Hi)(CR.XR.VM,  1,5.1  | 

CALL  GAUSS  I I ON .0 . 0 J I6M  .0 • .VI  I 
CAI  L uAUSSI  IFN.O. 0026.0.  ,V2 I 
('.At  L UAUSSI  IQN  . 0.  OJ  360 .0.  .V3  ) 

VMI  I I VP  I 1 ) A VI 
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YRI2  >=YR(2)FV2 
YR1 3) = VR<  3 I ♦ V 3 
CALL  (.MPROtC.X.Y  .3.9,  I > 

OO  80  I 3 1.3 
bj  yoi i ) *v«c  n -y ii  ) 

fUM  9TH  ORDER  MODEL  WITH  II  MACHINE  PAP AML  f F MS 


OU  70  I - I . I 
00  70  J = I . I I 
70  SYT<  J , 1 > = SY( l . J) 

CALL  OMl'MO  < S V T .WINV.SYrH.I  I .3.3) 

CALL  GMPROl SYTV.SV.SUM.I I .3.1 t > 

OO  79  1=1.11 

00  79  JM.II 

79  rscM<  i . j i * rsuMi  i . j) ♦ sum<  i , j) 

CALL  GMPRDISYTM.VO.YW.I  I ,3.1  > 

OO  89  1=1,11 

89  TYMl  I > = T YM ( I ) ♦ YM < I > 

Mol  CONTINUE 

CALL  M | M V I T SUM .1  1 .O.LMTX , MMT  X ) 

call  <; mproi  tsum,  tvm.delp.  1 1 . 1 1 . i » 

DC  AO  I 1=1.11 

PP  CjA  I N-ABSI  0£LP(  I > /PN<  I > ) -PL  I Ml  I 

IE  I PP  uA  I N.  (•  T . 0 . I DEL  P I 1 ) = PN(  I >•  (DELPI  I > / 

• AUSIOELPI I I ) I APL I MI T 
PNI  1 I =PN|  I MOELPI  I ) 

401  CONTINUE 

EM  9IH  ORDER  MOOEL  » I TH  II  MACHINE  PAPAME  TENS  END 

OO  7 1=1,11 

>’PArlQ(  I ) *PN  I I ) /P  TRUE  (I  > 

7 continue 

OO  9 1=1.11 

9 ESTMAT1  1ANI  T.I  » *PHAT10«I  1 
NNN  CONTINUE 

NrlEAP  =NREP  I TF  I 
MR  I TE(6.9) 

9 FORMA  T(  I Hi  * I OK  • *RA  T1  O • /// I 

MR  I TE  (6.6)  ( (ESTMAT  ( I.J),J=  1.111.1  = 1 ..NRFAP) 

6 FORMA  T<  .'0*.  1 IF  7.  3 I 

IF  (NPARAM.GT.3)  OELF  AC  =- 0 • 2 
Pf  AC  TR=PFAC  TRfOELFAC 
EBB  CONTINUE 
STOP 
t'NO 


SUBROUTINE  UMPRD ( A ,B .R.N.M.L ) 
DIMENSION  a < I I .8  ( I > ,R  ( I ) 

1 R = 0 
IK  =-M 

DO  10  K--I.L 
IK  = IK  *M 
OO  10  J - I . N 
IR=  1R  ♦ I 
J I =J-N 
18=  IK 
N<  IR ) =0. 

DO  10  I = I . M 
J I = J I *N 
lb= I8M 

1C  R(  IR  ) *R  ( IR  ) »A  ( JI  I *0  I I 8 ) 

U TURN 

* NO 
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TOIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
SROM  COPY  BURNISHED  TO  DOC 


MjltWlJ  l,  I I Nl_  HlNVU.N.D.L.M) 

UMI  NSIlIN  At  I ) ,L  ( t ) «MI  I ) 

0=1.0 
NK  = -N 

00  ttO  K - I .N 

UK  :NK  ♦ N 

1 ( K 1 = K 
M ( K > = K 
KA  -UK  »K 
HlbAsAlKKI 
Of.  20  J-K.N 
I l -U*  I J-  I > 

OO  2 0 l-K.N 
I J = I / ♦ I 

i:  III  AHSIIIIbAI-AHSIAI  I J)  I ) 15.20  ,20 

lb  Ult.A-At  I I) 

L I K ) = I 
M( K | = J 
.*C  CbNTINUI 
J:U(K  > 

If  l J—  K ) 3S.JS.2b 
25  K I *K — N 

OO  30  I = 1 , N 
K I = K I AN 
HOI-  O = — A I K I » 

J 1 =K  J — K ♦ J 
A(  K I I -At  J I ) 

JO  A | J | I iHtILD 
IS  I -Ml K ) 

If  I l-K  » 45.45.  3H 

!'•  jP  = N4  | I - | I 
DO  40  J = I , N 
JK  -NK  ♦ J 
JlsJP*  I 
MULO=-A(  JK  > 

A ( JK  J -A  I J I ) 

40  4 1 JI  > =MDUD 
4*.  IFIR1GA)  46.46.4U 
46  0=0.0 

RE  T URN 

4 H OO  55  I = 1 • N 

If  I I-K  ) 50.55,50 

5C  I K = NK  ♦ I 

Al  IK  ) =A  I IK  ) /I  -HI  v.4  ) 

55  CONTI  Nllr 

00  65  l-I.N 
I K -'  NK  ♦ I 
MUI.O  = A | I K ) 

1 J = I-N 

00  63  J = I , N 

1 J - I J *N 

If  I I -K  ) 60.65.60 

<>•(  II  I J-K  ) 62.65.62 

6 2 K .1  - I I * I * K 

A ( I J | *Hl!LO*A  IKJMAI  I JI 
6 5 CONTI  NUI 
K J =K -N 
0.1  / 5 J-I.N 

•J-K I KN 

II  I J-K  | / o,  15,  ro 

r 3 »(*  J)-AtKJ|/UIGA 

/«.  CoNflNU. 

O -l»4tt  |C,A 
AIKKI-I. 0/0  IGA 
'10  CONTI  NCt 
K =N 

I v ••  AsIK-ll 

If  IK  » ISO, 150. 105 
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2S0U  COPY  f UKNISHJSD  IX)  LDC 


I 0 b I -l  (K  ) 

IM  I -k  ) 120.120.  1 08 

I 0 a Jtl-N*  (K-  I ) 

JH  -N*  | I - I > 

UU  I t 0 J -I  .N 
JK  - JO*  J 
HOL  O- A I JK  ) 

J I - J H ♦ J 
A<  JK  ) = -A<  Jl  ) 

1 I 'J  A ( J l ) :H()LD 
I 2C  J-M<*> 

I*  ( J-K  > l 00,  I 00. 125 
I 2b  A 1 =K -N 

DO  I 1C  1=1  , N 
K I -M  »N 
H(jLD  = A(K  I ) 

J I -<  I -A  ♦ J 
A(KI)  --A  ( J!  » 

I 10  A I J|  > -MOLD 
GU  10  ICC 
ISO  WE  TURN 
E NO 


SLHRCJUTINt  oAUbS I I T.S.AO.VI 

A-0.0 

DO  SO  1=1.12 
call  randui i x , i y • v > 

I x s | v 
SC,  A 3 A A V 

W = TA-fc.  O ) A S ♦ A M 
H L I UR  N 
t NO 


SUb ROD  T I NL  RANOU 1 I X , 1 V . YFL ) 
1 V - | XA6SS.19 
IT  I I V > S , ft  . h 

5 lr:|»*2l4UHl(,.M| 

6 YT  L * I V 

VT  L=YFL«  .46S66I  iC-9 
RE  TURN 
ENO 


function  t cnt  t,  x.k  ,noe  t ) 

DIMENSION  F(IO).TIIO) 

DIMENSION  AR(S.S)  .A(S.S)  , U < S ) • US I S ) 
COMMON  AR.A.U.US 
FIR  1*0. 

IF (NOE T. tO. I » GO  TO  110 
DO  IOC  J = l • S 

10',  F IK  >=F(K  >AA(K,  J)AXI  J) 

GO  TO  I'iC 
110  DO  120  I 3 1 « S 
I 2 £ FTR)*F|A  )AAR(K,I)AX(|) 

ISO  IF (NOLI .EO. J)  GO  TO  300 
200  FlK)*H»  > ♦ L ( R ) 
v.U  TO  4<;c 

-ICO  FTK  >3F(F  >AUS|K) 

•CC  FCN=F(K) 

RE  TURN 
END 
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